THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING TRANSIENT HEAT TRANSFER TO A LAMINAR TUBE FLOW OF A BINGHAM PLASTIC WITH INTERNAL ENERGY GENERATION Carl D. Henning A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Mechanical Engineering 1965 July, 1965 IP-713

ACKNOWLEDGEMENTS The author wishes to express his appreciation to all those individuals who aided in the completion of this investigation. The following are especially acknowledged: Professor W. J. Yang, Chairman of Doctoral Committee, for his guidance throughout the course of the investigation. His criticism and suggestions were most helpful, and his encouragement aided greatly in the completion of the work. Professor J. A. Clark for his interest and pertinent advice. Professor H. Merte, R. B. Keller, and F. G. Hammitt, members of the Doctoral Committee, for their assistance, interest, and cooperation throughout the entire work. The Michigan Memorial Phoenix Project for providing financial aid. The Department of Mechanical Engineering for providing funds for the experimental investigation. The Computing Center, Dr. R. C. F. Bartels, Director, for providing computer time. The Industry Program of the College of Engineering for aid in the preparation of the manuscript. Finally, my wife for her sacrifices and understanding. ii

TABLE OF CONTENTS Page ACKNOWLEDGE MENTS.............................................. ii LIST OF TABLES................................................. LIST OF FIGURES.. a..*...........................vi NOMENCLATURE................................................... ix CHAPTER I INTRODUCTION................................ 1 II BINGHAM PLASTIC FLUIDS............................. 4 III ANALYSISa.. e.e........e.e....e e e.e *. e.. e 8 A. Step Responses to Inlet Fluid Temperature and Internal Heat Generation.......................... 13 B. Steady and Steady-Periodic Solutions to Internal Heat Generation........................... 18 C. Heat Transfer Computations**,* ****O,.............. 27 IV NUMERICAL RESULTS................................ 30 A. Steady Solution.*..................... 30 B. Internal Energy Generation Step Response.......... 35 C, Inlet Temperature Step Response................... 42 D. Steady-Periodic Solution to Internal Heat Generation*..................................... 46 V EXPERIMENTS............... *.... *................. 59 A. Experimental Apparatus............................ 59 1. Test section............................... 59 2. Fluid supply and heat exchanger..........*.... 62 5. Thermocouples................... 64 4. Heat generation device........................ 66 5. Instrumentation............................. 68 6. Bingham-plastic fluid....................... 69 B. Test Procedure.................................... 69 1. Velocity measurements......................... 70 2. Measurements and evaluation of physical properties of Bingham-plastic fluid........... 75 iii

TABLE OF CONTENTS CONT'D Page 3. Transient and steady state heat transfer 80 measurements...... o................... o o C. Results and Discussion....................o 81 VI: CONCLUSIONSo o o... o o.......... o......oo o o 92 APPENDICES................................................. 94 A DERIVATION OF TEE GOVERNING DIFFERENTIAL EQUATIONSo.. 95 B FINITE DIFFERENCE EQUATIONSo......................... 98 C ELECTRICAL ENTRANCE LENGTH ASSOCIATED WITH ELECTRODES e o............................................ 115 D CONFIDENCE LIMITSo.. o....o....................... 117 E COMPUTER PROGRAMSo................................ 121 BIBLIOGRAPHY o o o. o. o o o. o o o o o o o..... o o o. 134 iv

LIST OF TABLES Table Page I Internal Eenergy Generation Step Response Finite Difference Equations.......................... 16 II Steady-Periodic Finite Difference Equations........ 24 III Steady-Periodic Stability Requirements............ 26 IV Experimental Results for 25 per cent H2S04..... 82 V Experimental Results for 13 per cent Al 0 in H S04.. 83 VI Uncertainties in Variables........................ 118

LIST OF FIGURES Figure Page 1 Various Non-Newtonian Fluids....................... 5 2 Physical Model Showing Velocity Distribution and Elements in the Fluid and Tube Wall...............*. 9 3 Velocity Profiles for Various Values of the Bingham Plastic Constant'a'..........**...... 11 4 Radial and Axial Grid Used in Finite Difference Equations. *...................................... 15 5 Steady Fluid-Wall Temperature Difference versus Radius for Various Axial Distances.....*............ 32 6 Steady Centerline and Wall Temperatures versus Axial Distance as Function of a'...............33 7 Steady Fluid-Wall Temperature Difference versus Radius as Function of'a'........... 34 8 Step Response of Centerline Fluid Temperature versus Axial Distance.....................*........ 36 9 Step Response of Fluid and Wall Temperatures versus Radius...*...*......a....................... 37 10 Step Response of Fluid-Wall Interface Temperature as Function of'a'.........5.......... 39 11 Response of Fluid Temperature Gradient at FluidWall Interface Due to Step in Energy Generation.... 40 12 Step Response of Fluid-Wall Interface Temperature as Function of A and Fluid-Wall Heat Capacity Ratio M............................ 41 13 Centerline Fluid Temperature Response Due to Inlet Temperature Step........................... 43 14 Fluid and Wall Temperature Response Due to Inlet Temperature Step*................................. 44 15 Response of Fluid Temperature Gradient at FluidWall Interface Due to Inlet Temperature Response.,.. 45 vi

LIST OF FIG:RES CONTID Figure Page 16 Response of Fluid-Wall Interface Temperature Due to Inlet Temperature Stepo....................... 47 17 Amplitude-Ratio Response of Centerline Fluid Temperature as Function of a o o o o 48 18 Phase Lag Response of Centerline Fluid Temperature as Function of taT oO oO O O...... 49 19 Relationship between Sinusoidal Energy Generation and the Time at Which, a Fluid Particle Passes Location X = O O O O O O O O O O O O O O O O O 4 O O O O Q O O O O 52 20 Fluid Amplitude Ratio versus Axial Distance times Frequency for Various Radii O O O O.... o o o o o o o 54 21 Fluid Phase Angle Lag versus Axial Distance times Frequency for Various Radii o o o o.o.,......... 55 22 Fluid Amplitude Ratio versus Radius for Various Values of rX.......O...000000000000000OOOOOOOO...OOO 56 23 Fluid Phase Angle Lag versus Radius for Various Values of rx........... oo oooo..........oooooooooo^ 57 24 Periodic Temperature Gradient in Fluid at FluidWall Interface o.oooooo.ooooooooooooooooooooooo 58 25 Schematic Diagram of Experimental Equipment........ 60 26 Experimental Apparatus Showing Test Section, Knife Switches and Four-Channel Sanborn Recorder.. a 61 27 Detail of lest Section and Mixing Baffles........ 65 28 Friction Factor and Total-Static Pitot Tube Efficiency as Function of Reynolds Number o.......... 71 29 Radial Profiles of the Velocity Head in a Horizontal Tube o o OO o. o. O O O. O o o o o o 72 30 Relationship of Torque and Rotational Speed Obtained with Brookfield Rotational Viscometer for Various Bingham Plastics o 75 vii

LIST OF FIGURES CONT'D Figure Page 31 Flow Diagram for Various Bingham Plastics.......... 76 32 Rheolgi al Data for Slurries of Alumina and Water )......................... 77 33 Thermal Conductivity of Alumina Suspensions (5)..... 79 34 Step Response of Centerline Fluid Temperature....... 85 35 Steady Fluid and Wall Temperatures as a Function of Axial Distance....*........ a..,............... 86 36 Regimes of Free Forced and Mixed Convection for Flow through Horizontal Tubes (27)................. 88 37 Vertical Temperature Profiles in a Horizontal Tube.. 90 38 Vertical Temperature Profiles in a Horizontal Tube.. 91 39 Comparison of Numerical Results with Exact Solution by Michiyoshi, et al.(7)....................... 113 40 Test of Truncation Error Due to Axial Grid Dimension............................................... 114 41 Two-dimensional Infinite Channel with Source at Origin and Sink at + oo......................... 116 42 Computer Program A Flow Diagram, Step Response...... 122 43 Computer Program B Flow Diagram, Steady-Periodic.... 127 viii

NOMENCLATURE a Bingham plastic constant defined in Equation (3.2) A Coefficient of Sin rF in the periodic part of the dimensionless temperature in the fluid A Area, as defined locally in Appendix B; Aout, area through which heat flows out of an element; Ain, area through which heat flows into an element; Ain area of an element in the wall; Af, area of an element in the fluid; ft2. AR Amplitude ratio; ARf = +A7 + B2/E f; AR = =C' + D2/e dimensionl es s B Coefficient of Cos rF in the periodic part of the dimensionless temperature in the fluid C Coefficient of SinFG in the periodic part of the dimensionless temperature in the wall Cf Fluid specific heat, BTU/lbm-~F Cw Wall specific heat, BTU/lbm-~F D Coefficient of Cos rQ in the periodic part of the dimensionless temperature in the wall d Diameter of bob in Brookfield UL Adapter, cm D. Inside tube diameter, ft E Test section voltage, volts f Friction factor, hf/(L. u2 ), dimensionless g Local acceleration of gravity, ft/sec2 Gr Grashof number, p 9T4 T-T) L3/ L, dimensionless hf Fluid velocity head, ft of water I Test section current, amperes k Thermal conductivity; kf, fluid; kl, liquid; ks, solid; kw, wall; kwire, thermocouple wire ix

Le Hydrodynamic entrance length, ft 1 Length of bob in Brookfield UL Adapter, cm N Rotational speed, rev/min n" Slope of logarithmic plot of torque versus rotational speed, N, dimensionless p Pressure, lbf/ft2 Pr Prandtl number, JCf/kf, dimensionless q" Volumetric heat generation rate, BTU/ft3-hr q Mean volumetric heat generation rate, BTU/ft3-hr o ~q Surface heat flux, BTU/ft2-hr Q Fluid flow rate, ft /sec Qt" Dimensionless energy generation term defined following Equation (3.6) R A specific radius; R, radius of a specific point in the fluid; Ri, inside tube radius; Ro, outside tube radius; Rw, radius of a specific point in the wall; Rwire, radius of thermocouple wire, ft Re Reynolds number, puDi,dimensionless r Dimensionless radius, r*/Ri r* Radius, ft r Radius at which yielding takes place in a Bingham plastic, ft y s Ratio of bob to cylinder diameter of Brookfield UL Adapter S Specific gravity T Temperature; Tf, fluid temperature; Ti, initial temperature; To, entrance temperature; Tw, wall temperature; Tam, mixed mean temperature, OF t Time, sec u Dimensionless fluid velocity, u*/2u* m u*^ Fluid velocity, defined in Equation (3.3); ur*, average fluid velocity; u*ax maximum fluid velocity; ft/sec x

x Dimensionless axial distance x* Axial distance, ft y Dummy variable, dimensionless z Dummy variable, dimensionless o0< Fluid thermal diffusivity, ft /hr of/w Wall thermal diffusivity, ft2/hr le Dimensionless parameter, pfCf/pwCw Temperature coefficient of volume expansion, 1/ F r Dimensionless frequency, aiRi2/af Dimensionless wall thickness, -- 1 Ri ~C Dimensionless maximum periodic heat generation amplitude, defined in Equation (A-l) [ ~Volume fraction of solids in a suspension o Dimensionless time, t Cf/Ri A Dimensionless parameter, aw/af /1 Fluid viscosity lbm/ft-sec iU B Bingham plastic viscosity, lbm/ft-sec ^ Dimensionless temperature (T-Ti)/(To-Ti) 7 3. 14. Fluid density, lbm/ft5 Of Wall density, lbm/ft3 T ~Shear stress, T, yield shear stress of Bingham plastic; TwY sh"ar stress at wall, lbf/ft2' Phase lag of temperature with respect to the heat generation function; Of, phase lag of the periodic dimensionless temperature in the fluid; 0w' phase lag of the periodic dimensionless temperature in the wall; radians Y'e Potential function, ft2/hr xi

I fDimensionless temperature T-T)Kf -.. Ri2 yF. Dimensionless fluid temperature,t Dimensionless mixed mean temperature 3XwI Dimensionless wall temperature pf VPeriodic part of /f ltpW,, Periodic part of l 1Sf Steady part of If,ljS w Steady part of ~Jw (W Frequency, radians/sec xii

CHAPTER I INTRODUCTION Interest in heat transfer to slurries has been stimulated by the increasing emergence of non-Newtonian fluids as important raw materials and products in a large variety of industrial processes. With the advent of the homogeneous nuclear reactor employing the suspensions of uranium oxide or thorium oxide in heavy or light water as a reactor fuel, there has been research activity in steady heat transfer to slurries with internal heat generation. Slurries containing high concentrations of oxide tend to behave as Bingham plastics. Since the velocity at which the transition from laminar to turbulent flow occurs is rather high and erosion and corrosion become serious problems at high velocity in a typical slurry ), heat transfer to Bingham plastics in laminar flow is studied, despite the advantage of the high heat transfer coefficients realized in turbulent flow. (2) Poppendiek(2 has studied the fully-developed heat transfer both in laminar and turbulent flows of a Newtonian fluid through circular pipes with internal heat generation. Analytical solutions are given for the temperatures of liquid metals and ordinary fluids. Fluid temperature measurements, which are obtained in an experimental system with electrically generated internal heat sources, are compared with predicted values. Grigull has analyzed the fully-developed laminar heat transfer to non-Newtonian fluids flowing through a circular tube. Sparrow and Siegel( have made an analysis to determine heat transfer -1

-2characteristics for the laminar flow of a heat generating fluid in a circular tube with heat transfer. Internal heat generation is permitted to vary in an arbitrary manner both longitudinally along the tube and radially across the section, and longitudinal variations in wall heat transfer are present. Tachibana and Morishita(5) have experimentally investigated heat transfer in a slurry (consisting of distilled water and alumina) flowing through a circular tube. The experimental data was correlated by an equation similar to Dittus-Boelter's equation for Newtonian fluidso Schechter and Wissler) have extended the work of Sparrow and Siegel to include Bingham plastics, but only the case of an insulated wall has been treated. Michiyoshi, et al.(6,7) have further extended Schechter and Wissler's analysis to the case in which the tube wall is subjected to a heat flux. Solutions are given for both the entrance region and the fully-developed region. Ayers(8) has investigated the transient heat transfer from the wall with timewisely sinusoidal internal heat generation to the laminar flow of a Newtonian fluid in a circular tube. The asymptotic solutions are obtained by means of numerical computations. While steady heat transfer problems in slurries with internal heat generation have been fairly well treated, less is known about its transient behavior. This work investigates heat transfer to a laminar tube flow of a Bingham plastic with transient internal heat generation and inlet temperature. The nature of the disturbances is step and sinusoidal, with respect to time. An approximate solution is obtained by numerical methods using a set of finite difference equations derived

-3from the energy equations for the fluid and the wall. Both the fluid and the wall temperatures are functions of radial and axial distances as well as time. The experimental work includes the measurement of velocity and temperature distributions. The experimental results are in good agreement with theoretical predictions. The heat exchangers to which this analysis apply include the homogeneous nuclear reactor and a chemical reactor in which an exothermic reaction occurs.

CHAPTER II BINGHAM-PLASTIC FLUIDS Real fluids may be divided into two main classes: Newtonian and non-Newtonian fluids. According to Newton's law of viscosity a plot of shear-stress versus shear-rate for a given fluid is a straight line through the origin: i.e., the viscosity is a constant at a given temperature and pressure. Non-Newtonian fluids are those in which the viscosity at a given temperature and pressure is a function of the shearrateo This classification includes quite a few industrially important materials such as colloidal suspensions. Non-Newtonian fluids may be further classified according to the variation of the viscosity with the rate of shear. The shear-stressshear-rate diagram is presented in Figure 1 for various non-Newtonian fluids. The Bingham plastic, or ideal plastic, has a linear relationship between shear-stress and rate of shear, once it has been deformed. The relation may be written as A* T - (2.1) where 7T is the yield stress, the amount of shearing stress the Bingham plastic can withstand before deformation. The so called "real plastic"' has a constant viscosity at fairly high shear rate. The flow of a pseudoplastic may be described as shear-thinning and that of a dilatant as shear-thickening. -4-.

-5U, w. d u RATE OF SHEAR dru Figure 1. Various Non-Newtonian Fluids.

-6Colloidal slurries of solids in liquids have been found to be non-Newtonian by many investigators (9). Through surface-attractive forces the finely divided particles dispersed in the liquid are able to exert a strong influence on the mechanical properties of the material as a whole, when the particles are sufficiently small or sufficiently close together. In aqueous slurries of uranium oxide and thorium oxide this colloidal behavior gives rise to a definite yield stress below which the fluid will not deform (28) Thus, the Bingham plastic model has been used to describe its rheological behavior. Experiments (2) show that for an aqueous slurry of UO2 Ty = K, 7l e = ear and for Th02 y =3 u t = lhe e where 7 = volume fraction of solids in the suspension = viscosity of water K1, K2, K5, K4 = experimental constants which depend upon partical sizeo For a 1.4 micron particle size of UO or ThO2 K1 = 150 lb/ft K = 1.8 1 150 2 =2 K3 = 60 lb/ft3 K4 = 0.8 Sufficiently small solid particles will diffuse throughout the liquid due to their Brownian motion in a normal gravity field, and such dispersions are referred to as sols. For example, in aqueous solutions the ThO2 particles in a sol are less than 0.05 microns, while particles whose specific gravity is near unity may be as large as 0.5 microns 25)

-7When small particles coalesce to form loose, relatively independent clusters of particles, they are referred to as floes and may exhibit increased settling rates (although flocs of ThO2 have hindered settling rates 10-50 times those expected for unflocculated slurries(28)). Either suspensions of flocs or particles which are large enough to settle are referred to as slurries and often exhibit colloidal (25) properties (25) The thermodynamically stable condition of lyophobic sols is a flocculated state, but electrostatic forces produced by ions which collect on the surface of the particles cause mutual replusion, inhibiting the formation of flocs and reducing settling rates. Flocs also can acquire an electrostatic charge which results in hindered settling(28). This effect may be intensified by the addition of an electrolyte to the slurry as suggested by the reduced yield stresses in a ThO2 slurry due to the addition of oxalic acid(29). For a given liquid and solid phase the two most important variables which affect settling rates and fluid yield stress are concentration and particle size. Tests indicate that the settling rate decreases exponentially with increasing concentration(25) and yield stress is inversely proportional to the particle size(28). Either increased concentration or decreased particle size will result in both reduced settling rates and higher yield stresses.

CHAPTER III ANALYSIS The physical system to be studied consists of an insulated horizontal tube through which a Bingham-plastic flows steadily with internal heat generation (see Figure 2). Heat generated in the fluid is uniform spatially but may vary in either a step-wise or sinusoidal manner with respect to time. Temperatures in both the fluid and wall are considered to be functions of radial distance, axial distance, and time. Since the steady solution has already been obtained by Michiyoshi, (7) et al.(7) for both prescribed uniform wall heat flux and uniform wall temperature, and insulated outer tube wall is considered in the model. Due to the linearity in the differential equations of the system, the principle of super-position may be applied to obtain the solutions for the system with more complicated boundary conditions such as a specified wall heat flux. Since the Bingham plastic has a yield stress "T Ty the fullydeveloped laminar velocity distribution will exhibit a "plug" in the center in which the velocity will be a constant. If the radius of the plug is ry T_ - ( > r ry) (531) and _ _ = CL (532) Tw Then Equations (2.1) and (3.2) yield the velocity distribution. For r* 0~< R< a -8

WALL k \\ FLUID qrgd I o ~~~U(r*) -////// dr* VELOCITY f l~~~y I I I~~~~~qr*) r" PROFILE u(r*) x d x Figure 2. Physical Model Showing Velocity Distribution and Elements in the Fluid and Tube Wall.

-102 U, _ | -a ) (3 3a) IJC8 2Ga r < and for a< < 1 Ri Ac * - 2 C+_ R; (53.3b) __8 2 C K) The velocity profiles are illustrated graphically in Figure 3 for several values of'a'. Also, by defining the mean velocity to be U=i1 J (r*) 27T rdrK and applying Equation (3-3), one finds X - f LA C a - /- + 3 (354) In the formulation of the problem the following assumptions are made: a) The fluid is a Bingham plastic and flows in an incompressible, fully-developed, laminar, steady manner. This assumption is valid if an entrance length of the order Le > 0.035 D Re is allowed, (20) Re < 2300, and there is no time variation of the axial pressure gradient in the fluid. b) Axial conduction of heat in the fluid and wall are negligible. When the Peclet number (Be Pr), which represents the ratio of energy transport by enthalpy flux to that by heat conduction, is greater than 100 this assumption is valid. (24) c) The physical properties of the fluid and wall are independent of temperature. Fluid viscosity, however, is known to be quite sensitive

-111.0,-a =0.0 ya- 0.25 0.8-a=0.50 r 0.75 -a=.oI 00 0.4- 0.20.2 0.4 0.6 0.8 1.0 r Ri Figure 3. Velocity Profiles for Various Values of the Bingham Plastic Constant a.

-12to temperature. Thus, a large temperature gradient across the tube can lead to a considerable distortion of the velocity profile due to a radial viscosity gradient. d) The velocity and temperature are symmetric about the axis of the tube. Free convection can invalidate this assumption at low Reynolds number when there is a large radial temperature gradient. The buoyant forces associated with free convection can cause a secondary flow which increases the temperature at the top of the tube and decreases it at the bottom of the tube. This secondary flow also may cause a mixing of the fluid which raises the centerline temperature above the predicted value while lowering the average wall temperature. Any settling of solids suspended in the fluid will cause the velocity distribution to be distorted due to a viscosity gradient. The velocity will decrease at the bottom of the tube and increase at the top of the tube. With the assumptions, the application of the first law of thermodynamics to the system shown in Figure 2 yields the following partial differential equations for temperature in the fluid and wall. The velocity distribution u(r) is given by Equation (3-53) (Details of the derivation are given in Appendix A). Fluid (O < r < 1) J _ _ I Q (5_5) ^~-,L c Lka r F + + r + (5) Wall (1 < r< ) a -A(^ / 1r (5-6) >><9 - 7X 1 ap + j- at> )0

-130 for step disturbance in inlet fluid temperature where Q" = 1 for step disturbance in internal heat generation E (;1rB 0for sinusoidal disturbance in internal heat generation The boundary conditions to be satisfied are ( )J,) r 0 (initial) {o0iQ r) = o (entrance) a { X) rj = 0 (symmetry) a r ( X)i; 0 tv t( X 0i) 1) (fluid-wall interface) rI^|Aylj Y = 0 (adiabatic wall) A. Step-Response to Inlet Fluid Temperature and Internal Heat Generation The step response in fluid and wall temperatures is studied for two cases: Case 1. The inlet fluid temperature remains constant, and energy generation is suddenly started in the fluid. Case 2. The fluid has no energy generation, but the inlet fluid temperature is suddenly raised to a new level and maintained there. For a step disturbance in inlet fluid temperature the dimensionT T. less temperature symbol pf is replaced by f = T —- The equation for the fluid then becomes for the fluid then becomes -se * x'3 3X L r ar

-14The equation for the wall remains the same, but the entrance boundary condition is changed to.Ijo),G l = 1 These changes are so minor that it is unnecessary to treat Case 1 and Case 2 as two distinct problems. Instead, Case 1 is considered below, and important differences between Case 1 and Case 2 are noted. 1. Finite Difference Equations To solve the wall and fluid equations a radially and axially divided grid is constructed as shown in Figure 4o All derivatives in Equations (355) and (3.6) are approximated by finite differences involving the variables of surrounding points. The resulting equations are solved using numerical methods on the IBM 7090 digital computer. The finite difference approximations to Equations (3.5) and (3.6) appear in Table I. (Full details are given in Appendix B)0 Because the coefficients of the variables in the finite difference equations of Table I form a square 3(M + N + 1) tridiagonal matrix, where M and N are the fluid and wall divisions, respectively, the method of Gaussian Elimination can be used for their solution. Thus, an explicit method is used in the axial direction, while an implicit method is used in the radial and time directions0 This is implemented by solving for the variables at the first axial location with Gaussian Elimination. The axial dimension is incremented and the process continued until the desired axial distance is reached. 2. Computer Solution Computer program A, presented in Appendix E, is used to solve the finite difference equations of Table I on the IBM 7090 digital

Figure 4-K Radial and Ax~ial Grid Used in Finite Difference Equations.

TABLE I INTERNAL ENERGY GENERATION STEP RESPONSE FINITE DIFFERENCE EQUATIONS' Computer Eqn. Funct. Notation Location Difference Equation a. 4f(O,J) S(O) Centerline ( + VMAX - + 4N2) S(O) - 4N2 S(1) = + - VMAX U(o) + 1 AG L AG L b. lf(I,J) S(I) General Fluid S(I 1) + (VEL(I) + 2N2) S(I) — ) S( 1) +L( -N2(1 + 1) S(I + 1) = (I + P VEL(I) U(I) + 1 21 AG L c. (N,J) S(N) Fluid-Wall -(N - ) S(N-l) + ( + + (N- 1) + + 1) S(N) VW(O J) V(O) Interface 2 2 + ) - ( + 1) v(1) + ( + 1) F(N) B B 2 2N 2AQ BM N O d. 4w(K,J) V(K) General M Wall -A M -- MK (V(K-1) + (L + 2 2) V(K) 82 285(1+ -) j 62 L[ 28(1+- - ), ] e. Vw(M,J) V(M) Outside M M + 8 Wall A (+ M- 1) V(M-) + A( + V(M) = (AM+) W(M) * Dimensionless computer program nomenclature used in table is given in Appendix E Dimnsinles cn~nute p~g~a noen~~t~~e sedin a'ae i gienin Appendix E.

-17computer. The Program is written in the Michigan Algorithm Decoder Language. The Gaussian Elimination method used in the radial and time directions applies only to the solution of a tridiagonal matrix as would be produced by the set of equations below; where an, bn, and cn are coefficients of the variables yn' and dn is a constant. +I, + C Y d+ 2 Ad ^4SY + C3 23 = d2 (3.8) CC3 i; + $3 S C3 24 = d3 aM I1 #tIIYM = d The variables, yn, of the above system of equations can be found directly. Multiply the first equation by -a2/bl and add it to the second equation. The second equation is replaced by this resulting equation in which the coefficient of y1 is zero. Next, the third equation is replaced by the result of adding to it the second equation multiplied by -a3/b2. This manipulation eliminates the variable y2 from the third equation. The process is continued until every equation, except the last, contains only two unknowns. The last equation will contain only one unknown, making the value of YM directly available. By back substituting the value of YM into the M-l equation, whose only unknowns are YM and YM-1, the value of YM-1 can be found. Thus, by continuing the back substitution process the entire system of

-18equations is resolved. When made suitable for automatic machines, the algorithm for Gaussian Elimination is given by(26) M~~~~IA~~~~~~~ M ~(3o9) 44.^ ~i = M-l1 M-2,.. 1 where 7 and are determined by a, CC7/ A - cL, Ci = 2, 3,.-. M cL/ -- -~ _ ot c ----- i = 2, 53,. M This method is unconditionally stable if the equations are consistent. In practical applications this will always be the case, since each equation will contain a unique set of unknowns not repeated in any other equation. However, considerable roundoff error can accumulate because of the large number of arithmetic operations involved~ Therefore, it is necessary to check the final result. The transient solution should converge to the same steady-state solution already obtained by Michiyoshi (7) in closed form. If it converges to the correct steady-state value, the solution can be assumed to be accurate at all values of time, B. Steady and Steady-Periodic Solutions to Internal Heat Generation The term "steady-periodic" means that although the temperature varies periodically with time, the effects of an initial time condition

-19are no longer present. Thus, if a disturbance in the system varies periodically with respect to time, so the temperature will also vary periodically with a frequency equal to that of the disturbance. Because the governing equations and associated boundary conditions are linear, the solution may be divided into two portions: one a steady solution and the other a steady-periodic solution. The method of handling the steady-periodic solution is similar to that (8) used by Ayers ( Fluid |x,e, G ) = SF|X, r) r + Up4Xj)lr) ) (3-10) Wall W( X) r swx + ~ (Epw xa)r] (35.11) Substitution of Equations (3.10) and (3.11) into (3.8) and (3-9), respectively, and separating the steady and steady-periodic portions yields the following equations. Fluid YPS) _' Z: + 1 _ i. ) Wall C U s L - I'P a - o (3.14) r'- r r

-202/PWW A(+ +- Pw) ()315) The boundary conditions are kS;(o,)r) = z- #pC4oor) o t^f(x, o) _ _ P c ^^ (x e~ol o t'si(i) - SW( = xel - sw xo IR/P ) _ _ 0(Xe Ro/R)j O sinusoidal, q(l + Sin r ), it can be assumed that the fluid and wall temperatures will also be sinusoidal. Thus it is assumed ppoh [axFl + andxr)] c s, rte-n ) (5316) and iYnPs = [c (r) + D'x.r)]2 Sin( rc - b ) (t.17) where 1 -'A); = Tr (3 6)

-21By defining the amplitude ratio as the amplitude of the steady periodic temperature divided by the temperature resulting from an infinitely long period (F —O), the following definitions are derived. tA\~~R ^ —-78ZC~~;-~ ~(3.18) ARC w ~e (3-19) Then, the temperature response for the fluid and wall can be expressed as follows: = A |is + c ARI s (r0-F))e (3.20) K- % w t la ACRw S;( re - W) (3.21) It may be noted that the functions ARf, ARw, Of and..w are functions of x and r. 1. Finite difference equations The solution for the steady temperature distribution is obtained by applying finite differences to Equations (3.12) and (3.14) and solving the resulting equations on the IBM 7090 digital computer. Similarly, the steady-periodic temperature distribution is obtained. Equations (3.16) and (3.17) are applied to Equations (3.13) and (3.15), respectively, before they are written in finite difference form. For example, substituting Equation (3.16) into (3.13) results in

-22An cosrF - rsir + ~()(-A S g + -c CosFe] (3522) -j 3,A a %B p. J^B. I +I bAr e ar FY CoseS +S1.) r r Because Equation (3,22) is valid for all values of time, it can be divided into two equationso One equation corresponds to the sum of all the terms with Sin rQ coefficients, and the other corresponds to the sum of all the terms with Cos rQ coefficients. In this way the following equations are obtained. fr~A _ WA I DA _er + (r) aX = ^,-^+e Br + g3~23) D X a) r"- r r + Similarly the equations for the wall coefficients C and D can be obtained by combining Equations (3.15) and (3517) and separating the result into its SinFG and Cos Pr components. The process yields cr -f A(3DB ^ 3r + r r) +3o25)?cTr -_ *r + r dr ar - C<26) To solve the wall and fluid equations all derivatives are approximated by finite differences involving the variables of surrounding points. The finite difference equations for the fluid and wall are presented in Appendix B and the results given in Table IIo These finite

-23difference equations compose a set of simultaneous linear algebraic equations whose coefficients form a square 3(M + N + 1) matrix, where M and N are the wall and fluid grid divisions, respectively. An implicit method of solution is used in the radial direction, while an explicit method is employed in the axial direction. This permits use of an iterative technique by reducing the number of variables to be solved for at any one time to those present at one axial location. Starting at the first axial location, all variables are solved for by means of the Gauss-Seidel iterative method until convergence is achieved. Then the axial dimension is incremented and the variables solved for at the new axial location. This process is continued as far down the tube as desired. 2. Computer solution The computer program used to solve the set of equations in Table II is presented in Appendix E, and is based on one developed by Ayers (8). This program is written in the Michigan Algorithm Decoder Language for use on the IBM 7090 digital computer. The flow diagram also given in Appendix E schematically describes the order of computation used. Stability tests are automatically made to insure convergence for a particular set of data. If stability is not assured, the data is rejected and the next set of data is read. For the Gauss-Seidel iterative method, a sufficient condition for convergence requires that the absolute value of the element on the main-diagonal of the matrix, constructed by the coefficients of the variables, be greater than or equal

-24TABLE II STEADY-PERIODIC FINITE DIEREENCE EQUATIONS* Computer Eqn. Function Notation Location Difference Equation a. ^sf S(O) Centerline S(O) 4N2 S(1) + P/L VMAX U(O) + 1 P/L VMAX + 4N2 S(0 = AA 4N2 A(l) + P/L VMAX E(0) + 1 B(o) + e b. A(O,J) A(O) Centerline A(O) = N2 A(1) + P/L VMAX E(O) + B() +. B0J B0 Centerine B 0) 4N2 B(l) + P/L VMAX F(0) - F A(O) 4N + P/L VMAX c. B(O,J) B(O) Centerline B() B(1) P/L V F(O) - r A(O) 4N2 + P/L VMAX N2(l+ r) S(I+1) + N2(1-:) S(I-1) + r VEL(I) + 1 d.'('f(I,J) S(I) General S(I) = + / E Fluid 2N + P/L VEL(I) e. A(IJ) A(I) General N2(l+.) A(I+1) + N2(1- 1) A(I-1)+ VEL(I) E(I) + rB(I) + e Fluid A(I).-. -- 2N2 + P/L VEL(I) f. B(IJ) B(I) General N2(1+ 21 ) B(I+1) 21 1 ) + VEL(I) F(I) - rA(I) Fluid B(I), 21 21 -- 2N2 + P/L VEL(I) g. *sf(NJ) S(N) Fluid-Wall N(l+ ) S(I+1) + N(- ) S(I-1) + VEL() + 1 sw(OJ) V(O) Interface S(I) = V(O) = sw I 2N2 + P/L VEL(I) h. A(NJ) A(N) Fluid-Wall A() C(o) -(N - )(B(N)-B(N-1)) + ( + )(D()-D()) C(O,J) C(O) Interface A(N)= C(O) =-(N N N-) B 5 2 ( + 2) i. (NJ) (o) Ierfae B(N) FD(d) -A/B (M/s + 1/2) [C(l)-C(O)] + (N- 1/2) [A(N)-A(N-1)] - e/2N D(O,9J) D(O) Interface B(N) =D(O) r(5/2BM + 1/2N) J. isw(K,J) S(N) General Wall tsw 3 S(N) k. C(KJ) C(K) General A rM2 M 1 2A M2 A M2 M Wall C(K) -= I + 1 D(K+l) -2 D(K) + + S KS D(K- 1) r [5 25(1+ w)r 2 25(1+ - ) 1. D(KJ) D(K) General A M M 1 2 M Wall D(K) - = 5 25+ K) C(K+i) + 2 C(K) 2A M M A fg2 V- M C(K-l) j5 25(1+ ) m. C(MJ) C(M) Outer C) -A(M/5 + M- 0.5)(D(M) D(M-l)) wall C(M) Fs r M (1+6) n. D(MJ) D(M) Outer ) A(M/5 + M- 0.5)(C(M) - C(M-l)) Wall D(MC) O e r M (1+8) *Dimensionless computer program nomenclature used in table is given in Appendix E.

-25to the sum of the absolute values of all the other elements in the same row(26) Thus, for the i-th row: 3(M+-N -1) a-j > W |Qsi (35.27) j>1 % L The strict inequality must hold for at least one row. The stability requirements given in Table III can be determined directly from the equations in Table II. Inspection reveals that the equations corresponding to the steady solution are always stable. Inequalities d, e, and f are the most strict and are tested by the computer program. First, inequality d is tested at location N-1 where it is most severe. If it is not satisfied, the quantity P/L is increased until convergence is assured. An upper limit is put on P/L so that excessive computer time will not be required. Second, inequalities e and f are tested. If they are not satisfied, the data is rejected and the next set of data is read. Once stability is assured and iteration of the difference equations has begun, it is necessary to test for convergence after each iteration. Each variable is tested to see how rapidly it is changing with each successive iteration. For example, the variable S(I) is stored in location E(I) prior to each iteration. After each iteration of the wall and fluid difference equations, the following test is made. I l f- EIil < ALPHA.28) s(r) (j.28)

26TABLE III STEADY-PERIODIC STABILITY REQUIREMENTS* Inequality Equation Requirement a. Steady 4N Centerline - 2 + P 4N2 + L b. Periodic 4N2 + r Centerline 4N2 + VMAX L c. Steady 2N2 General 1 > Fluid 2N2 + VMAX L d. Periodic 2 + General 1 > — N+ Fluid 2N2 + VEL() e. Periodic 2(N - 1) + 2 (I + 1) Fluid-Wall 1 > + L Interface r(l + ) B 2M 2N f. Periodic 4AM2 General 1 > Wall rg. Periodic 2A(- + M - Outer 1 > Wall r (1 + b) * Dimensionless computer program nomenclature used in table is given in Appendix E.

-27The size of the parameter ALPHA controls the degree of convergence attained during the iteration process at each axial location. Errors resulting from incomplete convergence of the solution will influence the accuracy of the solution at the next axial location. Although the error may be small at the first axial location, the error will propagate in the axial direction and may become excessive. By varying the magnitude of ALPHA, it was found that a value of 0.001 was quite satisfactory. Decreasing ALPHA to 0.0001 resulted in a maximum change in the solution of 0.0025 per cent after 150 axial nodal points. When the convergence test similar to that given in Equation (3.28) is passed for all variables, the iteration process is ended. The program then continues to do other calculations before going on to the next axial location. The closer the initial guess is to the true value of the variables, and the further the inequalities (Table III) are from being equal, the faster the solution will converge. An upper limit of 100 is placed on the number of iterations, so that unnecessary computer time will not be used. C. Heat Transfer Computations The temperature gradient in the fluid or wall at the fluidwall interface is important to any heat transfer problem. For the steady periodic solution, the wall gradient is evaluated at j/6 time intervals according to the equation r N [(A(N )-A(t - iM r t ^(N) - B(N-Ico srFe (35.29) r i

-28In the step response solution both the fluid and wall gradients were evaluated according to Equations (3-30) and (3.31) respectively. r 1,r = N/(s|N) -s(N-1)J (5.50) ___ Al (v - v0o) (3531)'r r=1 - To achieve accurate values for the temperature gradients at the wall, a very fine grid was necessary. This is often the case when using a finite difference technique, because the fluid and wall temperatures change very rapidly near the fluid-wall interface. The mixed mean temperature as defined by Equation (3.32) is also evaluated. 2 TrJ r, (d r r T w - J C= *(K..r.U. (5332) 2 1 K i c, Tl(r*) dr* Using dimensionless parameters and expressing the integrals over the two velocity regimes of a Bingham plastic yields the following result for the steady-periodic solution. f(inf rdi4ffe r #c t hi2 b oIL r (/2lm(o 0a (3-33) J -iit)2 r1r tf+J I el-nc 2tr- r) rbor Using finite differences this becomes

-29z=-i w a Ll 1 - 1 2N ( ~l AA S~) Fs BeO| c<C6) rIN- Cad3 c|f-T f N ( 5|1) + A r-G t B ig cos re) 1- 6( 1/c j ( 81)cs rBJ/( /41,3 ) The equation for the mixed mean temperature, as it applies to the step response, can be obtained by dropping all terms in Equation (3.33) which have a Sin rF or Cos rP coefficient.

CHAPTER IV NUMERICAL RESULTS The IBM 7090 digital computer, located at the University of Michigan Computing Center, was used to run the computer programs given in Appendix E. In order to conserve computer time, most runs were made using the physical constants B A and S which applied to a 25 per cent sulfuric acid solution, with and without fine alumina particles, flowing through a polyvinylchloride tube of the size used in the experiment described in Chapter IV. The parameter'a' defined in Equation (3.2),which describes the velocity profile of a Bingham plastic, was varied over the range 0. - 1. to study the effect of a Bingham plastic on the temperature distribution. The numerical results are presented in four parts: steady, internal energy generation step response, inlet temperature step response, and steady periodic response. The truncation error inherent to the finite difference technique is discussed in Appendix B. A. Steady Solution The steady solution gives the temperature distribution in the fluid as a function of x, r, and a. Volumetrically uniform energy is generated in the fluid, and the outer tube wall is insulated. Because the temperature is time-steady and the outer wall is insulated, the tube wall has a uniform radial temperature equal to the temperature of the fluid-wall interface. -30

-31Figure 5 gives the steady temperature distribution in the fluid as a function of radial distance for various values of x; the Bingham plastic constant'a' equals 0.25. It can be seen that the radial temperature profiles are continuing to develop at x =.01. Michiyoshi(7) has shown that this development will continue until x =.30, after which both the centerline and wall temperatures will increase linearly with respect to x. Figure 6 gives the steady-wall and centerline temperatures as a function of axial distance for various values of the Bingham plastic constant'a'. Even for very small values of x the centerline temperature increases linearly with respect to x. However, the wall temperature still is increasing non-linearly at x =.01. Figure 7 illustrates how decreasing the value of'a' increases the centerline-to-wall temperature difference. This increase is due to the higher velocity present at the center of the tube when'a' is small. The effect of an increase in'a' is much more pronounced at large values of the parameter. For example, there is a considerably greater difference between the lines a = 0.50 and a = 0.75 than there is between the lines a = 0.0 and a = 0.25. The radius of the plug portion in the fluid velocity profile is proportional to'a', making the area of the plug proportional to a2. This squared relation accounts for the increased effect of'a' on the temperature profile as the parameter becomes larger.

-0.5X= 0.001 -I.0 -.5 - it -2.0t X- 0.005/ 1 > -2.5 - /^ 0a=0.25 -3.0 - X 0 0. 010 -3.5, Ri Figure 5. Steady Fluid-Wall Temperature Difference versus Radius for Various Axial Distances.

5 a=O. 04'..2 5_ 3_ WALL TEMPERATURE 03-" 1 ~- a=O. 50 2 - / —0a=0.25 CENTERLINE TEMPERATURE --- a=o.o 0 2 3 4 5 6 7 8 9 10 X* (RP)- x103 Ri e r Figure 6. Steady Centerline and Wall Temperatures versus Axial Distance as a Function of'a'.

a =0.95 o | -1.0 - = 0.75 0 -1.5 -2.0 M-a:=0.25/ "^ ~ _ ------- - ------ * - """7 " +~~X =0.005 =:o.o -2.5 -3.C 0I i _ i i I it 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ( r/Ri) Figure 7. Steady Fluid-Wall Temperature Difference versus Radius as Function of'a'.

-35B. Internal Energy Generation Step Response Program A, shown in Appendix E, was used to evaluate the transient temperature response of the fluid-tube system due to the sudden starting of internal energy generation in the fluid. The fluid and wall temperatures are functions of S,X,,,A, a, and r. At large values of time, the solutions obtained from the step response were found to agree with those obtained from the steady portion of the steady-periodic solution, and the analytical solution given by Michiyoshi (7). Figure 8 gives the transient centerline fluid temperature as a function of axial distance and time. At small values of time and large values of x, there is no x-direction temperature variation. This is because the enthalpy flux, which started at location x = 0 when 0 = 0, has not yet traveled sufficiently far down the tube. The characteristic time required for the enthalpy flux to reach a given location is x/u, The development with time of the radial temperature profiles in the fluid and tube wall can be seen in Figure 9. The equality of heat fluxes at the fluid-wall interface requires that the slopes of fluid and wall temperatures at the interface be of the same sign and is expressed as 3 t4' fA q-w' ar ~ 4 r ^(4-1) Inspection of the curves in Figure 9 reveals that the derivative of the fluid temperature, with respect to radius, experiences a very rapid change near the interface. This change of derivative appears

-565 _ 8=0.182 =2.4 =oo A=0.06 _ ( =0.0 tO 3- r 1.0 = 00.007 3-^-<( 2^ / /^/ =0.005 =0.004 =0.003 e=0.002 /"^^^i\ Q=0 O.001 0 2 3 4 5 RX (Rer) X 10 Figure 8. Step Response of Centerline Fluid Temperature versus Axial Distance.

3.0 —-li-iie 8__ =00C =0. 182 /P=2.4 2.5- A /=0.06 =0.30 2.5 - a =0.0 X =0.005 2.0 - 8=0.05-----------— FLUID —- - - - WALL — 11~ ~ ~ ~ ~ ~~~/s00.02 0.5 (r/Ri) Figure 9. Step Response of Fluid and Wall Temperatures versus Radius.

-38to be less drastic at smaller values of time. Evaluation of the radial temperature gradient of the fluid at the wall would require an extremely small grid size. Therefore, it was found more practical to evaluate the radial gradient of the wall temperature and relate it to the fluid gradient by means of Equation (41l)o An increase in the Bingham plastic constant,'a', yields a faster transient response of the fluid-tube system. Figure 10 indicates that the fluid-wall interface would reach steady state nearly twice as fast for plug flow (a = 1.0) as for parabolic flow (a = 0.0). This is due to the greater velocity present near the tube wall associated with an increase in la' for a constant mass flow rate. The radial fluid temperature gradient as a function of time is displayed in Figure 11. At zero time, it is equal to zero since the fluid and wall are the same initial temperature. As time increases, the temperature gradient goes through a maximum before decreasing again towards zero as steady state is approached. At infinite time, a zero fluid-wall temperature gradient is required because of the insulated outer tube wall. A few computer runs were made to determine the effect of the fluid and wall material properties upon the transient temperature response of the fluid-tube system. The increasing importance of the wall to fluid thermal diffusivity ratio,A, as the fluid-wall heat capacity ratio, M, is decreased is apparent in Figure 12. A high thermal diffusivity in the tube wall, large A, has a greater retarding effect on the fluid-wall interface temperature as the wall thickness is increased, M decreased. Heat conduction in the tube

1.0 0.4a =0. 95~ 0. 8 - =0..182 0.7 -a =0.50 07p /a=0.50-/ /\a =02.40 A=o.0 06.J -j 0.6- r=l. O: X / /X=0.005 4F 0.4- / // - 0.30.2 0. I I 2 3 4 5 6 7 8 F XeO2 Figure 10. Step Response of Fluid-Wall Interface Temperature as Function of'a.

-40-- a =0.0 a = 0.5 __- a = 0.5: 8 1.2 3 \ \ A= 0.12 X = 0.005 ~_2 \ X =0.00 3 2 4 6 8 10 12 14 16 18 20 8 x 103 Figure 11. Response of Fluid Temperature Gradient at Fluid-Wall Interface Due to Step in Energy Generation.

1.0 0.9 (P CpV )Fluid 0.8 M= (PCpV) Wall M=2II 0.7 - 0 0.6- - 0.4 0.3 - r i 32 0.2 - X = 0.005 0.2^/^ ^^ — A =0.06 0.1 I M= ---- A=30.0 0 0 I 2 3 4 5 6 7 8 0 xl02 Figure 12. Step Response of Fluid-Wall Interface Temperature as Function of A and Fluid-Wall Heat Capacity Ratio M.

-42wall results in a more uniform wall temperature response, retarding the fluid-wall interface response and increasing the outer-wall response. Greatly reducing the fluid volume or increasing the wall thickness, M approaches zero, would result in an infinitely long response time. C. Inlet Temperature Step Response A slight change in the data used in Computer Program A in Appendix E will enable the program, to be used for the temperature response of a fluid-tube system in which there is no internal energy generation, but whose fluid inlet temperature is suddenly raised to some new value and held there. If, in the input data, Q is put equal to 0.0 and INLET is put equal to 1.0, the inlet-temperature step response may be studied. The fluid and wall temperatures are functions of S X, X A,n a, and ro The temperature response resulting from a step increase in the inlet fluid temperature is displayed in Figures 13 and 14. As noted from Figure 14, for a = 0.5, the centerline response lags that for a = 00O, while the wall response leads. This difference results from the change in enthalpy flux due to the decreased centerline velocity and increased velocity near the wall as tat is increased. The radial fluid-temperature gradient evaluated at the tube wall is displayed in Figure 15. Initially the gradient is zero because the fluid-tube system is at a uniform temperature. When the enthalpy flux resulting from the increase in inlet temperature reaches a given axial location, the gradient increases, goes through a maximum, and

-431.0 0.9 \ N 0.8 0.7 - \< \ \ 0. 010o. \, 0.6 \ 0.5 ooo \ 0 1 =0.0052 4 - ^-x0.40 0.3 - =0.182z \=2.4 0.2 A =0o.o06'=0.001 — \ r = 0.0 I ~'-' --—'-a=0.5 I.a oi 0 I 2 3 4 5 X- (Re Pr)-x 103 Ri Figure 13. Centerline Fluid Temperature Response Due to Inlet Temperature Step.

0.9 ~ =0.015 -e=0.02 0.80.7 - 0.60=0.01-<Y \ \\ \ \\ \ o X ~~~80.04 8 —--- --- -8 \ 0.20 0.4~....... ~e,0.3-018 = o. -o /:2.4 A = 0.06 0.2 - X = 0.005 --- a = 0.5 K K 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.182 (r/Ri) Figure 14. Fluid and Wall Temperature Response Due to Inlet Temperature Step.

2.0 I I I i l l l l l I l l 1.8- X 0. 00 =0. 182 /=1.2 1.6- V^^A=o. 12 1.6'~\ ~:=-'-X =0.003 1.64 ~X =0.005 //.2 ~~ /*^ ^^I ---— ~ a =a =0.0 1.0 i T- 1.2- 0.2 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 & X103 Figure 15. Response of Fluid Temperature Gradient at Fluid-Wall Interface Due to Inlet Temperature Response.

-46decreases towards zero as steady state is approached. The peak value of the temperature gradient decreases with axial distance because some cooling of the fluid has taken place due to heat transfer from the fluid to the tube wall. An increase in the Bingham plastic constant raises the peak value of the fluid-temperature gradient. The greater heat flux due to the increased gradient is the factor which causes the more rapid temperature rise of the fluid-wall interface for larger values of ga?, as shown in Figure 16. D. Steady-Periodic Solution to Internal Heat Generation The steady-periodic solution gives the temperature distribution in the fluid and wall due to sinusoidal internal energy generation, qo c Sin Q@, in the fluid. Quantities considered are the amplitude ratio, AR =, and the phase angel, j. The steady-periodic solution is a function of g, X,,,,, a, and r. A resonance phenomenon is evident in Figures 17 and 18 which display the amplitude ratio and phase shift, respectively, versus axial distance times frequency. This resonance phenomenon is known to be characteristic of distributed parameter systems under distributed disturbances(12) As X increases, the amplitude ratio decreases; and the phase lag tends toward i:/2o Consideration of a simplified form of Equation (3-13), which assumes plug flow and neglects the effect of the tube wall, reveals some interesting results. The simplified equation becomes B - /- - r (4.2)

1.0 0.8 08.60X =0.003 0.6- / X=O.005 / / =o,=0.182 I / 0.4 / / -, I~-F.- // / / A=O. 12 - o I I I I r =1.0 0.2- / o=0.0 0 2 3 4 5 6 7 8 0 x102 Figure 16. Response of Fluid-Wall Interface Temperature Due to Inlet Temperature Step.

1.0 a=o 0.8- a0.25 = 0.182 a=1.0 \\\\a 0=0.50 0=2.4 0.6 \\\-a= 0.75 A =0.06 a 0=.95 0.4 - 0.2 0 0 2 4 6 8 10 12 14 16 18 20 rx Figure 17. Amplitude-Ratio Response of Centerline Fluid Temperature as Function of'a'.

2.60 -- 1 1 1 I — 1 a s 0.25 ___ - f r 2.40 - a 0.50 2.20 - a =.95 2.00 1.80 1.60 C,, z 1.40 1.20 0.80 - =0.182 = 1.0 0.60 -3 2.4 0.40 - A= 0.06 r 0.0 0.20 0.0 L 2 4 6 8 10 12 14 16 18 20 r x Figure 18. Phase Lag Response of Centerline Fluid Temperature as Function of'a'.

-50The boundary conditions are P (XO) = 0 (initial) (O, Q) = 0 (entrance) After solving Equation (4.2) and considering only that portion of the result which corresponds to the steady-periodic solution (Q > ), the following result is obtained. = -r -C.(12 rX-cos re + Si a SiPx ri e (4.5) Now for the case of r= 0 _ E z:X (4.4) ~r-o = 0 U The ratio of the steady-periodic temperature to the temperature when 0= can be written F ] _ l(i- Cos. rxCos re + Si 2PX Six r Ad T o2 rx - 21I-Cos 2P) [Fe -k T 1{- Cos 2 PrA (4 5).2 PX S;. 2 P, Equation (4.5) and Figures 17 and 18 show that for the case of plug flow (a = 1.0) there is a resonance phenomenon of period rX = I. As seen from Figures 17 and 18, decreasing the value of t'a increases the period of resonance until at a = 0.0, (Newtonian fluid) the period is about 2to The quantity rX used for the abscissa in Figures 17 and 18 was shown to be meaningful by the natural occurrance of rX in Equation (4.5). Also the plots of AR versus FX from numerical

-51solutions obtained for various values of r fell on the same line. The resonance phenomenon can best be explained physically by considering a single fluid particle as it travels through the tube. If the fluid particle passes the location x = 0, where heat generation is assumed to begin, with some inlet lag [r' (see Figure 19), a simplified equation can be obtained for an observer located on the particle by eliminating the space derivative of Equation (4.2). i = e Si(r e +re P (4.6) Integrating Equation (4.6) subject to the boundary condition Y (o) =0 yields #i = ( (i-Cosr CosrP' + sP oQ si Pe') The ratio of the steady periodic temperature to the temperature when 7= 0 can be written _ i-cosr) s 4r e 4, -i+ T-Co1gren (4.7) re[ (4.7) From Equation (4.7) it can be seen that the resonance phenomenon which appears in Figures 17 and 18 is due to the inlet phase lag shown in Figure 11. At [~ = 2t, 4t,... 2ng, the total energy generated in the fluid is equal to the energy which would be generated by steady state energy generation ( P= 0). The transient component of the particle temperature would be expected to be zero. However, the effect of the wall causes some transient component of the temperature to remain even at rG = 2i, 47,... 2nj. For plug flow the time Q can be replaced by the physical time x = 2x. Again it is shown that resonance occurs for plug flow with a period f = 2 or for plug flow with a period of rQ = 27t or PX =.

FLUID PARTICLE PASSES X= O HEAT GENERATION = E Sin(r) 4 \ -reo / l l - -E INLET LAG.....lI..... I I I I I 0 7T/2 7T 37T/2 27T 57T/2 37T rF Figure 19. Relationship between Sinusoidal Energy Generation and the Time at which a Fluid Particle Passes Location X = 0.0.

-53Figures 20 and 21 display the effect of radius upon the amplitude ratio and phase angle respectively. As the wall is approached (r ->1), the damping effect due to the wall increases. This causes the amplitude ratio to decrease and the phase angle to approach t/2. The amplitude ratio and phase-angle lag versus radius are given in Figures 22 and 23, respectively. A non-monotonic nature of the curves for larger values of OX is illustrated. Figure 24 presents the temperature gradient in the fluid at the fluid-wall interface as a function of time. The time-average of the wall gradient was evaluated numerically and found to be less than -9 10 9 A zero average wall gradient is necessary because the outer wall of the tube is insulated.

1.0 I r= 0 0.8 \ / r=0.5 8=0.182 \ \^ // ~=~~~~~E=1.0 -32.4 0.6 -r = 0.75 A=0.06 0.6 106 s- \ \\ r = 0.9 a =0.25 0.4- r=I.o 0.2 0 2 4 6 8 10 12 14 16 18 20 r x Figure 20. Fluid Amplitude Ratio versus Axial Distance Times Frequency for Various Radii.

I I I I I I I I 2.4 r=o r=0.5 2.0 1.6 o0lr 1.2 8 = 0.182',*F — gr 21 Fli Phas AgeLgeru~xaDitneims E= 1.0 Various Rad ii 1U.' -9- o.e- =2.4 A= o.o6 0.4 - = 0.25 0 2 4 6 8 10 12 14 16 18 20 r'X Figure 21. Fluid Phase Angle Lag versus Axial Distance times Frequency for Various Radii.

1.0 i I i I i 0.9 ~~~0.8- SL ~ = 0.182 6/ = 1.0 0.7 /3 —---------------------- -^p = = 2.4 A = 0.06 0.6 a = 0.25 r0.5 0.4 \ 0\ / —rx=4 _ 0.3 — ~~^^^-~^ / —-^x= \ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/Ri Figure 22. Fluid Amplitude Ratio versus Radius for Various Values of rX.

2.00 I I.80 1.60 - 1.40- rx=12.0 / iIA.o I |.00- I.X0 -r' x 0.40F- i. F- =2.4 0.20- Ao0.06 a=0. 25 0 0.1 0.2 0.3 0. 4 0.5 0.6 0.7 0.8 0.9 1.0 r/Ri Figure 23. Fluid Phase Angle Lag versus Radius for Various Values of [X.

2.0 1 8=0.182 r=8000 X =0.01 1.0 E=I.o oe o.oe _A=0.06 e' <a =0.0 -1.0 sinr Pf Or -2.0 -r /2 7r 37r/2 27r wt Figure 24. Periodic Temperature Gradient in Fluid at Fluid-Wall Interface.

CHAPTER V EXPERIMENTS An experimental program was conducted in order to determine the validity of the assumptions made in the analytical model, to ascertain under what conditions the analytical solution can be used, and to provide a confirmation of the numerical results. A. Experimental Apparatus The physical system to be experimentally investigated consists of an insulated horizontal tube through which a fluid, with and without solid suspension, flows steadily with electrically-produced internal heat generation. A schematic diagram of the system is shown in Figure 25. Temperatures were measured radially in the fluid and at the fluid-wall interface at the top, bottom, and sides of the horizontal test section. Mass flow rate, energy generation rate, and inlet and exit temperatures were simultaneously measured. The experiments were performed in the Heat Transfer Laboratory, Department of Mechanical Engineering, the University of Michigan. As illustrated in Figures 25 and 26, the component parts of the apparatus are; test section, fluid supply, heat exchanger, direct current power supply, and appropriate instrumentation. 1. Test section The test section is a 3/4" I.D., 7/8" O.D., polyvinylchloride tube, 1 foot in length. The fluid flowing through the tube is 25 per cent sulfuric acid, a good electrolyte, through which an -59

SANBORN RECORDER * —-CONSTANT- HEAD TANK |~..OVERO l POWER e 2 - S c D SUPPLY lZ. — _~Z POTENTIOMETER KNIFE -l - - _ XCHANGER |SWITCHE AN TEST MIXING S m INSULATED CALMING SECTION SECTION BAFFLES OVERFLOW HEAT LOWE EXCHANGER TANK PUMP Figure 25. Schematic Diagram of Experimental Equipment.

................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................ -...........................................................................................................................................v.......................................................................................................................................................................................................................................................................................... II...................................................................................... 1"...",:::::::::::::::::::::::::7::,::" 7:1................................................................................................... -............................................................................................................................................................................................................................................. I............I 11.1.............................................. I- -...................................................................................................... I'll. I 11- 1............ I..''..''............................................ I.-.............. 1.1- 1.................................................................................................................. I................................................................... I'll.,............. --.................................. Figure 26. Experimental Apparatus Showing Test Section, Knif e Switches, and Four Channel Sanborn Recorder.

-62electric current is passed from two platinum electrodes to produce internal heat generation. The polyvinylchloride tube was selected because of its corrosive resistance against the acid, besides being an electrical and thermal insulator. Preceding the test section was a 5 feet length of tubing encased in styrofoam insulation. This served the double purpose of establishing a fully-developed hydrodynamic condition as well as eliminating any influence on the inlet fluid temperature resulting from radial conduction of heat. Following the test section, mixing baffles were installed in the tube just upstream from the outlet fluid thermocouple as shown in Figure 27. This thermocouple read the mixed mean temperature of the fluid leaving the test section. The baffles are made of 0.1 inch thick polyethylene disks of negligible heat capacity with alternately spaced holes. 2. Fluid supply and heat exchanger A steady-state flow of test fluid was established by a constant-head tank located about 15 inches above the test section. The flow rate was measured using the discharge from the test section which was pumped continuously from a lower collection tank by a 1/15 horsepower stainless steel pump. Before returning to the constanthead tank, the fluid was passed through a heat exchanger consisting of a 50 foot coil of 3/8 I.D. polyvinylchloride tubing submerged in a constant-temperature water bath. Heat generated within the fluid in the test section was removed by the low temperature coolant in the heat exchanger.

THERMOCOUPLE WIRE ENCASED GAS VENT IN POLYETHYLENE TUBE PLATINUM STAINLESS STEEL ELECTRODE / TUBE C I I OPOLYETHYLENE MIXING Fgr 27. Dti o Tes Seto n gBAFFLES TEST) V^-^^ n^ ELECTRODE CASE CERAMIC INSULATOR Figure 27. Detail of Test Section -and Mixing Baffles

-643. Thermocouples The locations of the wall and fluid thermocouples, all made from 0.002" diameter copper and constantan wires, are shown in Figures 25 and 27. One fluid thermocouple was located upstream of the leading electrode, one immediately after the mixing baffle, and three at a distance of 7.8" from the leading electrode. The last three thermocouples were placed radially at the centerline and at distances halfway between the centerline and the inner surface, respectively. They were encased in 0.011" I.D. polyethylene tubes and placed in an 0.050" I.D. stainless steel tube, as shown in Figure 27. It was necessary to coat all thermocouples which extended into the fluid with General Electric 1201 Enamel. This coating protected the thermocouples from the corrosive effects of the sulfuric acid and eliminated any direct current voltage from being superimposed on the thermocouple emf. The final diameter of the coated thermocouple junction was estimated at.005 inch. Four wall thermocouples were installed at a distance of 7.8 inches from the leading electrode at the top, bottom, and both sides of the inner surface. Holes were made in the tube wall. Then the thermocouples were placed as near as possible to the fluid-wall interface and the holes filled with a solvent mixture of the wall material. Care was taken to electrically insulate the thermocouple wire from the fluid by an estimated.005 inch of the wall material. The inherent lag of the coated thermocouple made from 0.002" diameter copper constantan thermocouple wires was found to be less than 0.1 second by monitoring the transient thermocouple response with

-65a Sanborn recorder when energy generation was suddenly started in the test section while the fluid was stationary. By considering the thermal resistance of the air outside the tube, the tube wall, and the fluid inside the tube, the calculated heat loss from the tube was about 1 BTU/hr-ft2 for a 10 F temperature difference. This heat loss amounts to less than 0.1 per cent of the heat being generated in the tube. However, when considering the thermal resistance of the fluid and a 0.005 inch covering of wall material on the thermocouple, a 0.45~F error in temperature measurement results. Embedding a thermocouple wire of high thermal conductivity in a low conductivity material can cause an error resulting from increased heat conduction along the wire. J. V. Beck(30) considered the problem of a cylinder embedded in a semi-infinite solid perpendicular to a surface subjected to a heat flux. The equation for the disturbance is given below, where evaluation of the factor F was done numerically and presented in graphical form. (0) AT F (5.1) For a copper wire embedded a distance of two wire diameters below the surface of a semi-infinite polyethylene solid whose surface receives a heat flux of 1 BTU/hr-ft2, an upper limit on F may be obtained by extrapolating the graphs presented by Beck. Substitution of the appropriate values into Equation (5.1) results in a temperature disturbance of 0.25 F due to the presence of a wall thermocouple being embedded in the test section tube wall.

-664. Heat generation device Two platinum electrodes placed 11 inches apart at the sides of the tube were made of spirals of 36 gage platinum wire and enclosed in a polyethylene chamber as shown in Figure 27. A 3/8 inch hole was made in the tube wall next to the electrodes to allow electrical contact with the circulating fluid. The polyethylene chambers were located at a place where all gas generated at the electrodes due to electrolysis was caught and vented through a vertical tube. Thus, gas bubbles were prevented from entering the fluid flowing through the test section. The location of the electrodes at the side of the tube caused a non-uniform energy generation across the tube section, especially near the region of the electrodes. It is estimated in Appendix C that the electrical entrance length, defined as the one per cent variation of the potential across the channel cross-section, is about 1.3 diameters. This is not totally insignificant since the test section thermocouples were located roughly 10 diameters from the leading electrode. However, temperature measurements given in section C disclosed no consistent asymmetry in the fluid temperature resulting from the nonsymmetrical location of the electrodes. Because the electrodes were of finite size and the electrical potential was non-uniform near them, an effective length between the electrodes was taken as 10.5 inches. This length is subject to an uncertainty of + 0.5 inches. For an experimental fluid, a good electrolyte was necessary in order to electrically produce internal energy generation. A sulfuric

-67acid solution (25% by weight) was chosen because of its high (2) electrical conductivity(2). It is possible to eliminate electrolysis of the fluid at an electrode by using an alternating current power source of about (31) 400 cycles per second. However, when transient temperature measurements are to be made the resolution of the recording instrument is often reduced by the superimposing of an alternating voltage upon the thermocouple emf due to inductance between the power source and temperature-sensing instrumentation. Attempts to use both a 60 cycles per second alternating current power source and a full wave rectifier resulted in an a.c. pickup which was not eliminated by extensive electrical shielding. Finally, several 12 volt automobile batteries were used either in parallel or series for the power source. The voltage applied to the electrodes ranged between 24 and 48 volts. For the electrolysis of a strong acidic solution the reactions at the electrodes are given by(32) cathode reaction 2 H + e2 - 2 H anode reaction 4-oH — O2t + 2 20 + 4e then the net reaction is L H + a0-t 2 H4 $ 4t + 2 Ha~

-68Faraday's Law of Electrolysis states that the weight of a substance produced by a cathode or anode reaction in electrolysis is directly proportional to the quantity of electricity passed through the cell and is proportional to the equivalent weight of the substance. Using the heat of formation of water, Faraday's Law, and the definition of a faraday (magnitude of the charge of one mole of electrons) as being equal to 96.500 ampere-seconds, the energy absorbed at the electrodes due to electrolysis is found to be equal to 5.05 BTU/amp-hr. The resistance of the loop external to the test section was calculated to be at least 200 times that of the test-section, so that essentially no energy generation occurred externally. Considering electrolysis and resistance heating, the energy generation taking place in the test-section fluid is given by - 3Istum3n Et- ati(5.2) -71 R' z L 5. Instrumentation Instrumentation was required to measure and record fluid and wall temperatures, voltage and current, and fluid flow rate. Transient measurement of thermocouple emf and test section voltage were recorded on a four channel Sanborn 150 series recorder which had a maximum sensitivity of 10v/mm. The resolution of the line on the recording paper was estimated to be better than + 0.5 millimeters. A Leeds and Northrup model 8662 precision potentiometer was used to record the steady state thermocouple emfs. The resolution of the instrument was estimated to be 0.005 millivolts.

-69Volume flow rate was. determined by a calibrated 1100 cm3 vessel, accurate to 5 cm3, and a stop watch. The voltmeter and ammeter used to measure test-section energy generation were calibrated against Weston laboratory standard instruments. The resolution of the voltmeter was 0.2 volts and that of the ammeter was 0.02 amperes. 6. Bingham-plastic fluids As described in Chapter II, colloidal slurries of solids in liquids behave like Bingham-plastic fluids when the solid particles are sufficiently small or sufficiently close together. Slurries of titanium dioxide and Kaolin with 25 per cent sulfuric acid were found to be unsatisfactory because of rapid settling rates resulting from large particle size. Finally, aluminum oxide was adopted as a solid phase. After grinding the A1203 with water in a rotating ball mill for 10 days, the particle size was reduced to about 3 microns. On a microscope slide the particles formed flocs in concentrated areas, although independent particles were present as well. B. Test Procedure Major experimental programs conducted were velocity and heat transfer measurements for both Newtonian and Bingham plastic fluids, and the measurement of some Theological properties of Bingham-plastic fluids. Tube pressure drops and laminar velocity profiles were found with pitot-tube probes, while transient and steady-state temperatures were sensed by fine-wire thermocouples. Determination of viscosity and yield stress of A1203 slurries in H2S04 was done with a Brookfield rotational viscometer.

-701. Velocity measurements Attempts to measure velocity profiles and wall shearing stresses in slurries with pitot-tube probes were not very successful, Static pressure probes 0.10 I.D., spaced 50 diameters apart, were located in the tube wall. In addition, the radial position of a total pressure probe 0.030 in I.D. was controlled by a micrometer screw. A Chattock gage(33), precision manometer, was modified to allow measurement of differential water pressures by replacing the usual water phase of the gage with an oil of specific gravity equal to 1.2. This replacement increased the sensitivity of the instrument by a factor of 5, so that the resolution of the gage was 0.0013 inches of water. Prior to making tests with slurries, wall friction factor and pitot tube efficiency tests were conducted. The results are shown in Figure 28 where hf represents the velocity head. The measured friction factor agrees quite well with the usual Moody diagram while the totalstatic pitot tube probes are seen to be about 95 per cent efficient. Measurement of the velocity profile in distilled water yielded the expected results for a Newtonian fluid as given by Figure 29. The solid lines corresponding to the velocity headof a parabolic velocity profile were calculated from flow rate data. The tendency of the experimental points to lie above the calculated line could have resulted from the tendency of the observer to read the Chattock gage a little high. A slurry of A1203, 1,5 per cent by volume with a particle size of about 1 micron, was prepared and circulated in the test apparatus. After 20 hours of operation, a hard cake had formed on the bottom of the

0.4 -lI I I I 1|- I 4 0.2 - 64 2 Re 0.1 * _ ~-' —-----— ^\ v v v —: cu 0.08 - -0.8 CY NIEV'UE^ 0.05 0 FRICTION FACTOR f 0.5 _j _ — V PITOT TUBE EFFICIENCY,o b H ^ ~0.03 - -0.3 0.0 2 c rit -- -c0.2 0.01 ____ I I I 0.1I I 2 3 5 8 10 20 30 Re x1O02 Figure 28. Friction Factor and Total-Static Pitot Tube Efficiency as Function of Reynolds Number.

-721.0 0 H20, R, =1545 0.8 A 0 C H20, Re =1200 0. D 1.3 % A1203 - H20 Re 1180 0.6 - 0.4 -. CENTERLINE —-,, -' 0.2.4.1.2.3.4.5 VELOCITY HEAD, hf (in.H o) Figure 29. Radial Profiles of the Velocity Head in a Horizontal Tube.

-73tube which fractured when deformed (a report has been made of a 1/4 to 3/8 inch layer of ThO2 cake being built up on all parts of a 3 inch 200 gallon per minute circulating system(25)). Attempts to measure the velocity profile in the slurry were frustrated by the plugging of the pitot tubes with solid particles and time variant density changes of the fluid in the pitot tubes. The pitot tubes were down-flushed with water to clean out the Al203 particles and to remove any density gradient in the probes. However diffusion, small disturbances in the fluid flow rate, and radial movement of the total pressure probe caused particles to enter the pitot probes, cause density changes, and seriously disrupt the measurement of pressures. One velocity profile of uncertain accuracy is presented in Figure 29. The zero velocity at the bottom of the tube was known to be present by the inability of the pitot tube to penetrate to the tube wall. The velocity profile shifts upward, but does not exhibit the characteristic plug of a Bingham plastic and appears to be Newtonian in nature. 2. Measurements and evaluation of physical properties of Bingham-plastic fluids Increasing the concentration of the slurry to reduce settling and increase yield stress worsened the plugging of pitot tubes and inaccuracies of pressure measurements due to changing vertical density gradients in the probes. Finally a Brookfield LVT rotational viscometer with a UL adapter was used to Theologically define the properties of denser slurries. The UL adapter consisted of a rotating bob contained in an open cylinder with a ratio of cylinder to bob diameter, s, of 1.0985. By a method given by Metzner(19) the torque versus rotational

-74speed data in Figure 30 were converted into the conventional fluidflow curves shown in Figure 31. The shearing stress is expressed as _ 2 Tor e (535) The shearing rate was evaluated from the following equation: dr 4 1i/sNL d r 1 / 2/s [l -il+t K, ] (5 4) s2 I ( A l where K = 1 K S-l 2 652 and n" is the slope of the corresponding curve in Figure 30. The slurries used were relatively dense (13% by volume, Al203) and 95 per cent of the particles were less than 3 micronso Other investigations(4) of yield stresses and viscosities for slurries of A-203 in water are given in Figure 32. These results show rough agreement with the slurries of A1203 in 25 per cent H2S04 which were tested here. Experiments have shown that the evaluation of the thermal conductivity of a slurry by averaging the conductivities of the solids and liquid according to the volume fraction of each phase present gave predicted thermal conductivities which were erroneous by an order of magnitude (19) Orr and Dalla Valle(35) recommend the following equation for the evaluation of k of a suspension.

500 I I I I l 300 - 3'' 13 % AL203 IN H2S04 ) 100 - 40 - 0 10 % AL203 IN H2S04 DISTILLED H2 20 - 0.1 0.3 0.5 0.8 1.0 2 4 6 8 10 20 40 60 100 N [rev/min] Figure 30. Relationship of Torque and Rotational Speed Obtained with Brookfield Rotational Viscometer for Various Bingham Plastics.

4 3 13 % AL203 IN H2S04 E"~ IL/1B 4.06 cp i. 2. < 10 % ALtO3 IN H2 S 04 LB 3.26 p c ^^^ jy^^ DISTILLED H2O T e tL P 0.86cp 0 10 20 30 40 50 60 70 80 du' ( I/sec) dr* Figure 31. *Flow Diagram for Various Bingham Plastics.

/,~-) l'Xat pus'UT nurTV JO sa.mIJnTS IOJ'gI TeDISOTOGai'~ aJnS &'SaIlOS do NOII3OV8d 3Mnfl1OA g~'O ~'0 I'0 0 CA!' 3 I 09 09 OZHT./ A.JnlsTF OzH +4 ~O~Iv H08 ~~~~~~~~~~-, —------------------— 0

-78k, =, KsK K- - k I] 2 K K s I 77(N-Ksl where K1 and Ks represent the conductivities of the liquid and solid particles, respectively, and q represents the volume fraction of solids in the suspension. Whenever possible, experimentally obtained values for thermal conductivity should be used. Tachibana and Morishita5) measured the thermal conductivity of A1203 and H20 slurries. The results shown in Figure 33 were extended to slurries of A1203 and 25 per cent H2S04 as given by the broken line. For a 13 per cent by volume slurry of A1203 in H2S04, the thermal conductivity was taken as k = 0.52 BTU/hr-ft-F. The heat capacity of the slurry was taken as 0.1 W& AO + 0.81 Weh3t g (5 6) W ipktAg2o3 + WAIiLt I 2sok

1.0 0.9 X 0.8 0.7 0.6 A203 - H20, EXPERIMENTAL 0.5 o O 0.4 AL203 + H2S04' 0.3 - U 0.2 0.1 0 0.0 0.05 0.10 0.15 0.20 0.25 VOLUME FRACTION OF SOLIDS, r Figure 33. Thermal Conductivity of Alumina Suspensions.(5)

-803. Transient and steady-state heat transfer measurements The transient response of the fluid and tube temperatures was studied when subjected to a step change in the internal energy generation. A 25 per cent by weight H2SO2 solution with and without fine suspended alumina particles, 95 per cent less than 3 microns in diameter, was circulated through the test loop. Without the suspended alumina, the fluid corresponded to the limiting case of a Bingham plastic (a = 0.0), a Newtonian fluid. For 13 per cent by volume of Al203 suspended in the H2SO4 solution the Bingham plastic constant,'a', was calculated from rheological and flow rate data. Equation (3.4) can be put in the following form: ^t 4 C'm _ 0 - 40- 3 1. -4Q+3 ^(57) __Uo 12 CL Since the factor a' cannot be solved for explicitly, a trial and error procedure was necessary using the rheological data from Figure 31 and the flow rate data corresponding to a given test condition. Test section voltages ranged between 24 and 48 volts, while the Reynolds number was varied from 368 to 1333 for the H2S04 and from 417 to 685 for the A1203 - H2SO4 slurry. A typical value of internal energy generation, 80,000 BTU/ft3hr, resulted in a bulk fluid temperature rise of about 2~F.

-81All tests were conducted with the inlet fluid temperature at room temperature to reduce any heat transfer to the test section. The inlet temperature was adjusted by varying the heat exchanger coolant temperature and waiting whatever time necessary for thermal equilibrium to be established. Three flow measurements were made before each test. Once the desired flow rate was achieved, all thermocouple emfs were recorded with a precision potentiometer. The test section power was then turned on and transient centerline and wall temperatures recorded on the Sanborn recorder. Once steady state had been reached, the temperatures were again recorded with a precision potentiometer and three additional flow measurements made. C. Results and Discussion Experimental results obtained for a 25 per cent by weight H2SO4 electrolyte are presented in Table IV. These results correspond to a Newtonian fluid or Bingham plastic with constant a = 0.0. Similar results for a slurry of 13 per cent by volume of Al203 in H2SO4 are presented in Table V. The Bingham plastic constant, ta, ranged between 0.30 and 0.39 for the slurry data. While the internal energy generation calculated from the rise in bulk fluid temperature agreed to within 5 per cent of that calculated from test-section voltage and current for the H2S04, a disagreement was noted in the slurry data. Measured bulk fluid temperatures would indicate a greater energy generation than that given by the test-section electrical power measurements. Since this is physically impossible, the disagreement is attributed to the

-82TABLE IV EXPERIMENTAL RESULTS FOR 25% H2S04 Temperature (oF) To(inlet) 71.5 71.5 71.5 71.5 71.0 71.2 71 top-TO 4.44 7-77 5-55 5.69 6.09 5.16 6.89 Tenter-To 7.55 1.33 0.49 0.62 0.80 0.40 0.755 Tbottom-TO 1.38 2.76 2.13 2.00 2.98 2.13 2.62 T*ide To 1.51 3-33 2.27 2.09 2.40 1.91 2.93 Tside-TO 1.51 3.20 2.58 1.96 2.89 2.31 3-02 Texit-To 1.56 2.80 1.33 1.11 1.78 0.977 2.04 Volts 23.4 33.9 33.0 34.0 34.0 34.0 34.0 Amperes 1.24 1.95 1.95 1.96 1.99 1.97 1.98 Re 368 471 963 1190 790 1333 634 Pr 10.8 10.8 10.8 18 00.8 00.8. 10.8 X x 103 5.27 4.12 2.02 1.63 2.44 1.46 3.06 qtemp(BTU/hr-ft3) 34,400 79,400 77,000 79,300 83,600 78,200 77,500 qelc(BTU/hr-ft3) 35,140 81,730 79,460 82,400 83,660 82,820 83,240 x 102 steady aver.wall 2.19 1.78 1.34 1.24 1.47 1.19 1.58 steady centrln 0.678 0.513 0195 0.237 0.02 0.152 0.286 GrPr x 10 0.4 0.8 0. 08 0.8 0.8 0.8.8 * Wall temperature on same side that electrodes are located.

-83TABLE V EXPERIMENTAL RESULTS FOR 13% A1203 in H2S04 Temperature (~F) To (inlet) 78.4 78.4 78.6 78.6 78.7 top-To 2.22 2.66 5.11 1.55 1.60 Tabove center-TO 0.222 0 355 0.666 0.222 0.311 Tcenter -To 0.400 0.444 0.89 0.311 0.800 Tbelow center-TO 1.55 5.50 16.0 3-77 10.09 Tbottom-To 4.66 10.2 24.9 5.11 15.40 Tide -To 1.33 333 8.56 2.00 5.55 Tside-TO 0.89 2.99 5.04 2.22 4.12 Texit-To 0.579 0.890 2.04 0.489 0.756 Volts 24.2 36.2 48.3 24.3 36.2 Amperes 0.98 1.49 2.00 1.0 1.50 Re 417 540 533 530 685 Pr 7.22 7.22 7.22 7.22 7.22 X' 103 7.00 5.40 5.46 5.58 4.25 temp (BTU/hr-ft3) 38,600 76,600 174,000 41,500 82,600 qelec (BTU/hr-ft3) 28,784 66,867 124,467 29,881 67,320 ilf ~ x 102 I steady aver. wall 5.00 4.32 5.33 5.29 5.84 ~ steady centrln x 102 0.752 0.359 0.387 0.563 0.643 a 0.39 0.34 0.34 0.34 0.30 * Wall temperature on same side that electrodes are located.

-84filling of the cavities between the mixing baffles with solid particles resulting in inefficient mixing of the fluid. For each test the transient response of the centerline fluid temperature was recorded as shown in Figure 34. For 25 per cent H2S04 the experimental measurements give a faster response than predicted at the initial state, and then fall below as the system approaches its steady state. However, the discrepancy between the theory and experiment is within an acceptable range. For the A1203 suspension the measured response falls below the theoretical prediction but approaches the predicted steady state at large values of time. This deviation is thought to be due to a slower centerline fluid velocity than expected, resulting from settling of the solid phase or adhering of solid particles to the thermocouple probe. The 95 per cent confidence limits obtained by the method outlined in Appendix D are presented in Figures 34 and 35. The steady-state temperatures for the fluid and wall are shown in Figure 35. For 25 per cent H2S04, a Newtonian fluid, the experimental data of the centerline fluid temperature agree fairly well with the theoretical curve. However, a vertical assymmetry is noted in the wall temperature. This deviation is thought to result from free convection which was neglected in the analysis. The buoyant forces due to temperature differences in the fluid could cause a secondary flow which results in a radial mixing of the fluid, raising the centerline fluid temperature and lowering the average wall temperature. This conclusion is supported by Metais and Eckert's study(27) of forced, free, and mixed convection regimes for horizontal tubes.

4.8 I |I| 4.44.0 0O H2S04,X=0.002, aoO.O 3.6- V A203 + H2S4,X=0.00425, a =0.3 3.2 2.8X=0.00425- 2.4 —,., —, v.^0^^^ ^-X= 0.002 0.8 95% CONFIDENCE LIMITS — " X lo3 Figure 34. Step Response of Centerline Fluid Temperature.

-8620 I I 0 H2S04, a=0.0 VB VB 10 v A1203 + H2S04 VB 8 _0.3< a<0.4 VB VB 6 SUBSCRIPTS A- AVERAGE WALL VA V B- BOTTOM WALL A 4- C- CENTERLINE FLUID 0 A V T S- SIDE WALL VS TS. 3- T- TOP WALL - 0 2 - OT T VT - WALL OA ^~0 >-^OA~~~T 0 0 6 OB S S T B I O 0 S GB OB B B 0.695 % CONFIDENCE LIMITS OC'C 0.4 - vc CENTERLINE / VC 0.3-,,,, i 0.2 -_ _ oc. -- a=0.0 0.1 1 2 3 4 6 8 10 X (RePr) xl103 Ri Figure 35. Steady Fluid and Wall Temperatures as a Function of Axial Distance.

-87Their findings for the case of constant wall temperature are shown in Figure 36, where the Grashof number, Gr, is based on the tube diameter and the difference between wall and fluid bulk temperatures. These results can be compared to the case under study here: internal energy generation in the fluid with insulated tube wall. As the quantity X* (RePr)-L is decreased, the average wall temperature in Figure 35 Ri is seen to approach the predicted value. The quantity GrPr I associated with each test is nearly constant and equal to 8 x 103. Since X- (Pr)1 Ri is nearly constant, as Reynolds number is increased and GrPr D held L constant the effect of free convection is diminished. This would correspond to approaching the forced convection region in Figure 36 from the mixed convection region along the line GrPr = 8 x 103 This apparent agreement between Metais and Eckert's results and those presented here might have been anticipated because of the similarity of fluid temperature profiles resulting from either a heated tube wall or internal energy generation in the fluid and an insulated tube wall. For 13 per cent by volume alumina in 25 per cent H2S04, the centerline temperature agrees with the theoretical prediction. However, the bottom tube wall temperature is significantly higher than expected, causing the measured average wall temperature to fall above the predicted value. This disagreement was due to settling of the solid particles to the bottom of the tube. A non-symmetric vertical velocity gradient was established which greatly affected the temperature profiles. The denser, slower moving fluid at the bottom of the tube lost less heat due to the enthalpy flux, and became warmer.

-88_ i _-_ _ __ FORCED CONVECTION TUfBULENT FLOW |Id ==_;_ Nu=.a116[ +(-)] (Re2-125) Pr'() _ _- HAUSEN 10' - TRANSIT1O0 LAMINAR-TURB LENT i;>r-. — L — 1 1 1 1s, ^ ^^ I _ I I t I ITURBULENT FLOW ReD se - Nu=4.69Re P"'" (()~'" Re I.s3 l l l [. METAIS v I FORCED CONVECTION LAMINAR FLOW Nu1-.86 PeF/()" / MIXED CONVECTION,) SIEDERaTATE / | LAMINAR FLOW SIEDER a TATE t |Nu1.75. ) Gz +G0,008.O3(GrPrr] Pt OLIVER -— 1FREE CONVECTION HOLDEN - KERN o METAIS A PETUCHOV ---- x SHERWOOD 10 Io ido' id' id o' lo' io' i' o" Gr Pr DFigure 36. Regimes of Free, Forced, and Mixed Convection for Flow through Horizontal Tubes.

-89Further illustration of the vertical temperature profile can be seen in Figures 37 and 38. One curve in Figure 37 corresponds to the Newtonian fluid, H2SO4, and exhibits the increased upper tube wall temperature thought to be due to free convection. All the other curves for the A1203 slurry exhibit an increased lower wall temperature which is attributed to settling of the solid phase. A time dependent effect is disclosed by Figure 37, resulting from the continued settling of the solid phase. The test conditions corresponding to the triangles and hexagons are nearly the same. However, the asymmetry of the temperature profile is becoming more pronounced as slurry-circulating-time increases. The tests as listed in Table V were separated by approximately one-half hour intervals from left to right. Thus, the test corresponding to the hexagons in Figure 37 took place about 2 hours later than that corresponding to the triangles. From the vertical temperature profiles it appears that the solid phase is slowly continuing to settle and increasing the fluid temperature assymmetry. Examples have been cited where the solid phase of a slurry has continued to settle until the flow was entirely choked off(25). Even in turbulent flow, a laminar sublayer exists in which particles can settle(20). The settled particles may act as the effective tube wall, so that settling can continue until an equilibrium thickness is reached, or the tube is completely blocked. The rigidity of a settled bed appears to be inversely proportional to the particle size. Therefore, slurries utilizing large particles may form a vertical density gradient, but the thick mud at the bottom of the tube will continue to move rather than forming a hard cake and eventually plugging the tube.

1.0-y-""'0 H2S04, q"= 79460,X = 2.02 x 103 0.5 - to 3T~ 0. > 0'^ VA1203+H2S04,q" 28784,X=7.0 x103 0 A1203+H2S04,q" 29881,X=5.58x163 SUBSCRIPT'S' REFERS TO SIDE WALL TEMPERATURE 0 SS W s CENTERLINE | 5 0 R V K 1.0 I I 0 2 3 4 5 6 TEMPERATURE (OF) Figure 37. Vertical Temperature Profiles in a Horizontal Tube.

1.0 // V A1203 H2SO4, q"=124467, X=5.46x10'3 0 A1203+ H2S04, q=66867, X=5.40x10'3 0.5 0 A1203+ H2S04, q"=67320,X=4.25x10-3 ) II SUBSCRIPT'S' REFERS TO SIDE 5a * WALL TEMPERATURE C) U) S S $ -J Z \ s^ CENTERLINE | 9 5 0.5 - 1.0 L 0 5 10 15 20 25 30 TEMPERATURE (OF) Figure 38. Vertical Temperature Profiles in a Horizontal Tube.

CHAPTER VI CONCLUSIONS 1. The largest steady centerline-to-wall temperature difference exists for a Newtonian fluid, (a = 0.0). As the Bingham plastic constant'at is increased, the centerline-to-wall temperature difference decreases until at a = 1.0 (plug flow) it is equal to zero. This is due to the uniform velocity and insulated outer tube wall. 2. The response time necessary for the fluid-tube system to reach equilibrium is decreased by an increase in the Bingham plastic constant?a'. A resonance phenomenon is observed in the steady-periodic temperature response. The period of resonance equals 2t for a Newtonian fluid (a = 0,0) and equals t for plug flow (a - 1.0). The heat capacity of the wall dampens the temperature fluctuations of the fluid near the wall. 3. The wall heat capacity has a large effect on the transient temperature response. This effect is one of increasing the length of time required to reach steady state and is intensified as the ratio of wall to fluid heat capacity is increased. 4. The regimes of forced, mixed, and free convection proposed (27) by Metais and Eckert (27) for a horizontal tube with constant wall temperature, appear to be valid for the case of internal energy generation in a Newtonian fluid with an insulated tube wall. 5. Settling of suspended particles in a slurry distorts the velocity profile and produces a higher than expected temperature -92

-93at the bottom of a horizontal tube when energy is being generated in the fluid. Experimentally, the centerline fluid temperature agrees with the predicted value, but the average wall temperature is higher than expected,

APPENDICES -94

APPENDIX A DERIVATION OF THE GOVERNING DIFFERENTIAL EQUATIONS Under the assumption listed in Chapter III, the first law of thermodynamics can be written for the model showin in Figure 2. Fluid X t ax" ~,i 3T-~r Tr HE. S x +AI1X1 2 r 2 f - (A.1) t r/^ rp*2 r Dp* acf Wall aTw E DTTWY 1 4w + (A.2) t c: r r Dr*z For the three types of problems considered, steady-periodic, internal energy generation step response, and inlet temperature step response, the term q" in Equation (A.1) takes the following forms energy generation step: q" = q inlet temperature step: q" = 0 t > 0 steady periodic: q" = qo e sin 6t The boundary conditions for the general problem are T( X% 0, r) - T7( X O, rr) =- T (initial) Tk(Ot, ) ]= To (entrance) Tif(XtO) = o (symmetry) r 7i;(X t, )) - Tw(X t R')?(fluid-wall interface) aTV r( r,)? R (adiabatic wall) -95

-96By defining the following dimensionless variables l =- (T s-r)/:for disturbances due to internal heat o R;I generation r- T _= -T-:for step disturbance in inlet fluid temperature oL T- | r _____ _ (A.3) r RC'.2 at R; RI{'0 2 LRJ Equations(A.1) and (A.2) may be written in dimensionless form as Fluid! XX = B + { + i' (A.4) -x r z r %r Wall JQ~'^=f r + (A.5) ~ r r r (A 5 where Q" = 1:for step disturbance in internal heat generation Q" = C SinPe:for sinusoidal disturbance in internal heat generation and Q= A 0"= 0:for step disturbance in inlet fluid temperature. The boundary conditions are j(x,o,r) = 4)WxO, ) o J 0:for disturbances due to internal heat generation 1:for step disturbance in inlet fluid temperature only.

II 03 7 OII I c, Q"o Q) X:^ 1 -- — "I n-.- (-^ ri

APPENDIX B FINITE DIFFERENCE EQUATIONS 1. Notation As shown in Figure 4, the physical model in Figure 2 is divided radially and axially to form a grid, and a system of notation used by Ayers(8) adopted. The center of each element is referred to as a nodal point whose temperature is the average temperature of the element. Any nodal point in the fluid can be described by using a double subscript notation (I, J), where I is the radial counter and J is the axial counter. For a wall nodal point the radial subscript is K, so a point in the wall is referred to by the subscript (K, J). Radially, there are N + 1 fluid elements and M wall elements. It follows that the radial grid dimension in the fluid is 1/N, while in the wall the radial grid dimension is Ro/R; - 1 _ The radius of any fluid nodal point is given as I/N, (I = 0, 1, 2,... N), K5 and the radius of any wall nodal point is 1 + M,(K = 0, 1, 2,... N). Axially, there are P divisions of the tube length L, so the axial grid dimension is L/P. The axial distance is JL/P, which is measured from the point where energy generation begins, X = 0. 2. Finite Difference Approximation of Derivatives If a function u = u(x, y) is assumed to have enough partial derivatives, the value of the function at the two points (x, y) and -98

( A (x + ax, y) can be evaluated by the Taylor's series. L, + AX L + 2 (B.1) The above equation can be solved for the first derivative a =t + o (x) (B.2) xjr - ^x 4 ow (*2A 12 X Equation (B.2) is known as the forward difference approximation to the first derivative. The truncation error involved is of the order Ax. A more accurate approximation to the first derivative can be had by writing Equation (B.1) for the point u. Th s+ (B Then, subtracting Equation (B.3) from (B.1) and rearranging the result yields )(,t, * (af - ^d + ~ o{Ji)X (B.4) Equation (B.4) is the central difference approximation to the first derivative. It has a truncation error of order of magnitude x. The central difference approximation to the second derivative of u(x, y) can be obtained by adding Equations (B.1) and (B.3) and eliminating unwanted derivatives.

-100aZ2 " _I. A l, IL +* o IXI (B. 5) Whenever possible, the central difference approximation to derivatives should be used because of the smaller truncation error involved. 3. Steady Finite Difference Equations The partial differential equation for the sJteady temperature distribution in the fluid may be obtained from Appendix A as I $r +1 (B.6) <jr)~'- - -— r r + Note that the wall need not be considered, because at steady state it has a uniform temperature equal to the fluid-wall interface temperature. When writing the finite difference approximations to Equation (B.6), the axial counter need not be used since the resulting equations are only solved at one axial location at a time. This makes the designation of the nodal points more simple. Fluid Wall T5frj) = $r S ) t4ywr)j = v () 95(|r+1v ) = SIz+.I qs w K+I) = V K+1 Y8S l r -ilJ )= S Ir-1) TWI K-1 JI}= V IK-v ) SsF I( a-1) = - K|

Using the above notation, the finite difference approximation for the terms in Equation (B.6) can be written _S S - |ULIP (B.7) d/ Lp {'L/P I a> f = 1 / sl^ - s | P br ~ T /NI 2/N I (B.8) a t f s| | + sir-1) - 2 s( (B'r 2 i//N v (B.9) Note that a central difference approximation is used radially, but axially a forward difference approximation is used. This is made necessary by the explicit method used in the axial direction. Substitution of Equations (B.7), (B.8), and (B.9) into (B.6) gives the finite difference equation for a general point in the fluid. Fluid s(i) = N2l'- ) sIt*)+N/(A- -)sfr-i-)+ V~LI]) U(r) + (B. 0) 2NL + P VEL Il A special equation must be written for the centerline and fluidwall interface elements, as special conditions exist at such limiting nodal points. At the centerline, the term in Equation (B.6) sf-r is not defined at r =. By l'Hospital's rule it is found that

-102Lirr I _F _ _ __s r-o r r - Prz Therefore, the governing equation for the fluid at the centerline of the tube is ~rn aX _ -l ^^..11) Crux ~W>Xf - 2 3 ^ +~ (Bl) Using finite differences, Equation (B.11) becomes s4o) 4 / iV s Ls I VMA U J) 1 (B.12) P VMAX + 4 Nz The fluid-wall interface element, whose fluid subscript is (N, J), has a thickness equal to one-half of a general fluid element. The governing equation becomes o0 - KG r | + %~ v (B.13) In terms of dimensionless parameters, the preceeding equation becomes nR %=: rL R4 (B.14) Ri 3r rti R; where Re 1 S = 1 R. - 2, R- 2A

-105Application of finite differences then yields the following equation s \N-/2 + (B.15) 4. Steady-Periodic Finite Difference Equations The partial differential equations governing the steady-periodic temperatures in the fluid and wall have been found to be Fluid 1a o~U atirI'-i a'Kitf + _ a3/pft I d;- [s (B.16) Wall s'w A a ar 4 - (B.17) -r~ Let the periodic temperatures be written as follows pF - A(x,r) r S re + B(X,r)Cos Ce (B.l8) %^ -- C(X/r) Sin P0 i D(Lr)Cos ^ (B.19)

-1 4Substituting Equation (B,18) into (B.16) gives the following result for the fluid. Ar Cos [') - l(t rS,-,t- (- (l)/ -- ^ CrS o; 2-A D An / ^ - -'-P4r + c ", Fo +o, r.' C r " + E iB r20 (B.20) Because Equation (B.20) is an identity, it can be separated into its Sin rQ and Cos rF components. A B - - +2 + (B.21) DA yA ( ) -Br + r) X = - r t r F +- (B,22) Again the notation can be simplified. A|rT) Jr= A) B|ZIJ) = BII) Ar+1, J) - A(Ir+) B(|i+TJ - B( r +1 A(T-l,7]. A i-T] B(r- T) 1= 3 IB|-I) A IT T-1i = E(r) B(r,;-() = _ ( ) Employing finite differences with the above notation results in the following two equations for the fluid.

-105A() = ir'2r1 B(rl + (I- ) -) - r L - r CL BEJ V Z LF r(r (B,23).2N2 L V vE4I) L r rB B|I| = l Ai T + - Al ^- I) A f- ^ ~i N 4- L VEL II rJ + P VL I (B. 24) By the same procedure given in Equations (B.18) through (B.24), the following result can be obtained for the wall. r 2Y CIK)-+ +'g a 2(1+ L K- K (B. 25) IM _ I 12AM (B.26) The equation for the centerline element takes on a special form because a term in Equation (B.16) is changed,

-106L,' j _ zf. 2 f r 0 r ar 3 rApplying Equations (B.18) and (B.19) and finite differences to the modified form of Equation (B.16) yields the following result for the centerline element. A~o 0 ~ - F 1BiO 4-, j VMo — F,o (3.27) rr L ao(0) -A(l( r V ) j0 + V M (0)- [ 4 (B.28) Bo= r Al +. AO V- LrP 0 (:B.28) The fluid-wall interface element is a composite of both fluid and wall materials. The energy equation for this element can be written IPw v~c, -9 v~,,:.'rJT' - 4, 1^ ~c~ ^ v;)| 3-t = ^, A, -tA- t Y, cE s,,,,,,t (B.~) where VW = 2t R.w 2M - Kw, J| 7ro r As = rn R i +i Aot = 2r R(1+2 ) Now = - t Ai^ =.?7rR;(i-ZT-1

-107A combination of dimensionless parameters, the above quantities, Equations (B.18) and (B.19), and finite differences, results in Equation (B.29) becoming the following two equations.'AN | = { t21 l 1tst1 1\/V 11 + {g, {~l-{0|(B.30) ^)._ IN-.-^(S(^ -^^2l)M l B(AI) A/ z= N (B-1 + 24A 2N The element at the outside of the wall (M, J) has a thickness of one-half that of a general wall element. The energy equation is found to be WCWVw |it = 1 AN (B.32) where AIN 2 TrR e(I, ) Vw= 2 7 M AI' TT~(t —2?

-108Using the above quantities, dimensionless parameters Equation (B.19), and finite differences, Equation (B.32) is transformed into the following two equations. Ci'/ = A S M' -"- ~ (B,33) D(M|) - F -,X, ( ci -c - (B.34) 5. Step Response Finite Difference Equations Both the inlet temperature step response and internal energy step response are so similar that only the latter case will be considered. The only difference separating the two cases is the energy generation term. in the fluid equation which is not present for the first case. The energy equations for the fluid and wall were given in Chapter II. They are repeated here. -beT f lr) - t r a1f i (B.35) I ~ ~ r? dr\~C

-109Notation used in the finite difference equations can be simplified by eliminating the axial subscript and the superscript relating to time increments. This is possible because the variables are solved for at one axial location and for one time increment during any phase of the solution. Fluid Wall X+-1 1M+1 |I, (&| r, T) =w i K,J)9= V (KJ 5(ri*) v(l) u"'f [r-1 - S| K-1, V =/ K-1I I (r,-1) = UTJ = w(k) )'IJ|= F(l) Central finite difference approximations to the derivatives are used in the radial direction, but forward differences are necessary in the axial and time directions due to the method of solution. The terms in Equation (B,35) can be written as follows'%,sl.,)- F'FrJ P _ sLr-i) *SX ~ L/P

-110i ^4. | s _ I_ -sl-) |_ r r r/N 2/N',.3 I' "' I / 1V. - Substitution of the above quantities into Equation (B.35) gives -A/l. l' 11 + ( st.- ( -j V) ( -N 1-2-1js-.)s -- 4 L VELI U(F (. 37) The equation for a general point in the wall can similarly be found. Equation (B.36), changed to apply to a centerline element, takes the following form.'w~e-L + AX^' 2- (B. 39 In terms of finite differences, the centerline equation becomes

I++ v VMAX F4N2js(O) - 4/\z(1) = 0)+ -VMAX U(O) + 1 (B. 40) The composite fluid-wall interface equation appears as 2 M 2MJ ~ 2NW1 ril A, A/,,t 5 I bv| + (B,41) $+ IN,cg 2M) r I2N If it is noted that S(N) V(O), Equation (B.41) takes the following form when written as a finite difference equation. -(N- 2) S(N-1+(1 I # #)# 4 + (l)4N-2)s(N) -A & }2 2- v F2N (B,.42) The energy equation for the half element at the outside of the wall can be written __ 1 A / / =/ - g- IV, (B.43)!ll/ -l A1 + ~ ~.... (B. 43) ~ ~ ~~

-112Applying finite differences results in the following equation for the element located at the outside of the wall. -A(l+M-4 2)( V2M-l| V( M +A, ifM-l2l) VWM = Sl l (|Mj (B.44) 2 AOAI 6. Truncation Error Because the finite difference approximation to a derivative is 2 found from a truncated Taylor series, an error of the order (Ax) is introduced. This truncation error is directly dependent upon the grid size used in the numerical solution. Too few radial or axial divisions in the fluid or tube will cause a considerable deviation between the true solution and the numerical result. However, too many divisions will require more computations and increase the computer time. Figure 39 shows the comparison between the exact solution by (7) Michiyoshi, et al.' and the numerical solution with 20 radial grid divisions. Figure 40 shows the effect of axial grid size. It can be seen that convergence to the correct solution is quite rapid as P/L is increased. Evaluation of temperature gradients near the fluid-wall interface is difficult because of the large temperature changes which occur there. In the step response solution as many as 100 fluid and wall divisions were used to accurately evaluate the temperature gradients.

0=0.5 ^^ -0.01- X =0.009OSH I " MICHIYOSHI I " N=20 3o -0.03 -0.004 -00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 FoRi Figure 39. Comparison of Numerical Results with Exact Solution by Michiyoshi, _t al.(7

10 8 9- 8=0.182 x ARf =1.0 CJ 8- /32.4 7- A=o.06 I ~ r=6ooo 00 1 $N ~6 a -0.25 C\M \X =0.005 "a: 5 - r =1.0 4 sC 2\ 1 ^ 4~~~~~~~~~~~~~~~~~~ w~~~~~~~~~~~~~~~~~~~ DL 3 >^ 2 400 1200 2000 2800 (P/L) Figure 40. Test of Truncation Error Due to Axial Grid Dimension.

APPENDIX C ELECTRICAL ENTRANCEI LENGTH ASSOCIATED WITH ELECTRODES A non-uniform energy generation across the tube section in the region of the test-section electrodes results from their nonsymmetrical location at the side of the tube. An estimate of the entrance length required before the electrical potential is uniform across the tube can be obtained by considering a two-dimensional infinite channel of width Di which has a source at the origin and a sink at + co (see Figure 41). The potential for this configuration can be obtained by adding the potential which causes a uniform flow of gM/2Di along the positive real axis to the potential resulting from an infinite number of equi-distant sources of strength 2CM separated by a distance 2Di along the imagninary axis(56). The resulting potential, fe, is given by - h 4 ( cs- Co)- (C.1) e 2 Di 2 Di Numerical substitution reveals that at Z = 153D. there is only a one per cent variation of the potential across the channel section. -115

y/ S OC -. ST TS Diink - + 00 SOURCE OF STRENGTH 27TM Figure 41. Two-dimensional Infinite Channel with Source at Origin and Sink at +.

APPENDIX D CONFIDENCE LIMITS By a method outlined by Kline and McClintock37), an uncertainty interval or confidence limit can be placed on an experimental result which is a function of several variables that possess an uncertainty themselves. The uncertainty interval is specified as that interval about the mean into which a certain percentage of the values of a variable or result would have fallen if the experiment had been repeated an infinite number of times. When the result R is a function of n independent variables, then the uncertainty interval is given by where' V = indep n t + a (D.1) where: V = independent variables n WR = uncertainty interval of the result Wn = uncertainty interval of the independent variable Table VI contains typical mean values of important variables with an uncertainty of which the author is 95 per cent confident. These values will be used to determine the confidence limits for typical results in order to illustrate the method employed and the factors considered. Energy generation in the test-section fluid is seen from Equation (5.2) to be a function of test-section voltage, amperage, radius,. and length between the electrodes. The uncertainty interval -117

-118TABLE VI UNCERTAINTIES IN VARIABLES Description Mean Value Uncertainty (+) r* 0.5371 in. 0.002 in. Length between electrodes, L 10.5 in. 0.5 in. x* 7.8 in. 0.3 in. Test-Section Voltage 36 volts 0.2 volts Test-Section Current 2 amps 0.02 amps Position of fluid thermocouples 0.050 in. Absolute thermocouple reading 75~F 1~F Wall thermocouple error due to heat loss to ambient 10~F 0.47~F Wall thermocouple temperature disturbance 10~F 0.25~F Resolution of Sanborn Recorder and Potentiometer 0.2~F q" 88,915 BTU 4489 BTU hr-ft3 hr-ft3 Wall temperature 100F 0.84~F Centerline fluid temperature 0.5~F 0.284~F

-119of this result is given by W is 3 E I'~E) + WI + (B RI, SL WL (D 2) By taking the appropriate derivatives and using the values given in Table IV, the uncertainty in energy generation is found to be W //= 4480 BTU/hr-ft The mean value of energy generation is 88,915 BTU/hr-ft3. Thus, the uncertainty is about 5 per cent of the total. The centerline fluid temperature uncertainty is a function of the resolution of the recording instruments and radial and axial location of the thermocouple [ ~T 2 T2 I?T 21 WT [I DR WiR) + W rx + DXJ W (D-3) L)R Assuming the mean value of the centerline temperature rise to be 0.50F, the uncertainty may be evaluated T [( ( 2 (, S 2 (o, i) (0o,2) 4-. - - 0_ _,. 2 +.2\) W',T - 0,7/ 1 7 (D.4) This uncertainty equals 56.8 per cent of the total centerline fluid temperature. The uncertainty in the wall temperature measurement is a function of the resolution of the recorder, axial location, heat loss from the tube wall, and the disturbance due to the presence of the thermocouple itself. For a 10~F wall temperature rise, the uncertainty interval is 0.84~F, which amounts to 8.4 per cent of the total.

-120The uncertainty interval or confidence limit on the result can be expressed in a dimensionless form _ T krk (D.5) 2~1.i cr In~ i W I ( WR- K,]7 (D.6) \1/ T\ R I For the dimensionless centerline temperature olzi4= 4t SM ^+l 0 2j _2 0 ~$7/ (D.7) and for the dimensionless wall temperature ___ oP [ 44) 2 (48 A2(7) zl,./ S (D.8) The 95 per cent confidence limit placed on the dimensionless centerline and wall temperatures would be 57.1 and 9.85 per cent, respectively, of the total, based on the typical mean values considered.

APPENDIX E COMPUTER PROGRAMS The Computer Programs A and B were used to solve the finite difference equations listed in Tables I and III for the step response and steady-periodic temperature distribution, respectively. The steadyperiodic computer program and nomenclature are based on those developed (8) by Ayers () A flow diagram preceeds each program to schematically show the order of computations. -121

-122START _ READ DATA PRINT DATA INITIALIZE VARIABLES COMPUTE CONSTANTS MATRIX COEFFICIENT MATRIX INCREMENT TIME PRINT IRESULTS F INCREME ( J> P --- X T STEADY STATE F _ _- REACHED X -0 Figure 42. Computer Program A Flow Diagram, Step Response.

COMPUTER PROGRAM. NOMENCIATURE STEP RESPONSE A(I) Coefficient of first independent variable in Equation I ALPHA Convergence limit B(I) Coefficient of second independent variable in Equation I BETA pf Cf/Pw Cw C(I) Coefficient of third independent variable in Equation I CMPARE WGRAD(A/B) COUNT Time step counter DELT Time increment DELTA Dimensionless wall thickness, 6 D(I) Constant terms in Equation I FGRAD Temperature gradient in fluid at fluid-wall interface F(I J) Storage matrix for rf(I, J) FREQ Time frequency of solution printout I Radial counter in the fluid, (I - 0, 1, 2...N) INLET Inlet fluid temperature J Axial counter, (J = 0, 1, 2...P) K Radial counter in wall, (K = 0, 1, 2...M) L Length of grid network in the axial direction LAM A M Number of radial elements in the wall MIEDT Mixed mean temperature of fluid N Number of radial elements in the fluid -123

-124P Number of axial elements in the fluid and wall Q. Internal energy generation switch, (Q = 1.0,. there.is energy generation; Q = 0.0, there is no energy generation) RATIO Bingham plastic constant, a S(I) Fluid temperature, tf T Time, 0 U(I) S(I, J-l) V(I) Wall temperature, VEL(I) Fluid velocity, VEL(I) = u(r) VMAX Maximum fluid velocity,'vMAX = uax/2u WGRAD Temperature gradient in wall at fluid-wall interface W(IJ) Storage matrix for $w(K, J), x = (x*/Ri)(RePPr) X Axial distance

-125COMPUTER PROGRAM A - STEP RESPONSE DIMENSLON A(200),B(200),C(200),D(.200),S(200),EPS(200),GAMMA(200), *001 I F(6000,AC),W(4000,AQ),V(100), VEL(100),U(100) *01 VECTOR VALUES AC- 2,20,20 *o2 VECTOR VALUES AD- 2.20,20 *003 INTEGER COUNTFREQ,N,MIK,J *004 PRINT GOMMENT S2HEAT TRANSFER TO BINGHAM PLASTIC IN LAMINAR *005 I FLOW WITH INTERNAL HEAT GENERATION *** STEP RESPONSE ***S *005 START READ DATA *006 PRIOT~Of- MMENT S1 NEW DATA$ *007 PRINT RESULTS N,M,RATIOFREQ,DELT,BELTLA ELTAL,PALPHAQ, INLET *008 WHENEVER (P-N.G.6000).OR.(P*M.G.4000).OR.N.G.100.OR.M.G.100 *009 PRINT COMMENT $DIMENSION OF ARRAY EXCEEDED. DATA REJECTED$ *010 01 END OF CONDITIONAL *011 01 VMAX -(3.*(1.-RATIO).P.2.)/(RATIO.P.4.-4.*RATIO+3.) *312 WHENEVER RATIO.E.I. *013 THROUGH LOOPAt FOR I-0,1,I.E.N *014 01 LOOPA VEL(I)-0.5 *015 01 01 TRANSFER TO ENTRYA *016 01 END OF CONDITIONAL *017 01 THROUGH LCOPB, FOR I-O0,lI.E.N *018 IF-I *019 01 NF-N *023 01 VEL(I)-3.i*(.-2.*RATIO+2.eRATIOeIF/NF-IF.IF/NF/NF)/ *021 01 1 RATIO.P.4.-4.*RATIO+3.) *021 LOOPB WHENEVER IF/NF.LE.RATIO, VEL(I)-VMAX *022 01 ENTRYA CA - NON *023 CB a P/L *024 CC a N-0.5.025 CD a 0.5/OELT*(DELTBETBETA/M+1./N) 026 CE - LAM/BETA*(M/BETA+0.5) *027 CF - LAMNM*M/OELTA/DELTA *028 CG a 0.5*LAM*M/DELTA *029 CH -DELTA/M *030 Cl a LAMo(M/DELTA+M-0.5) *031 CJ - O.5OELTA*(1.+DELTA)/DELT/M *032 PRINT COMMENT $THE CONSTANTS GENERATED BY THIS DATA ARES *033 PRINT RESULTS CA,CB,CCCDCECFCG,GCH,CICJVEL())...VEL(N-1) *034 COUNT * 0 035 T a DELT *036 THROUGH LCOPD, FOR J-Ot,1,JG.P *037 THROUGH LCOPC, FOR O-0,1,IG.N *038 01 LOOPC F(I,J) * O. *039 32 THROUGH LCOPOD FOR K-1,1,K.G.M *040 01 LOOPO W{KJ) ~*. *041 02 B(O) - 1./DELT+VMAX*CB+4.*CA *042 C(0) a -4.*CA *043 THROUGH LCOPE, FOR It1,1,I.G.N-1 *044 All) ~ -CA*(l.-0.5/1) *045 31 1B() I 1./DELT+VEL(I)*CB+2.*CA *046 01 LOOPE CM(I) -CA(1.+0.5/1) 047 31 A(N) a -CC *048 B(N) a CD+CC+CE *049 C(N) - -CE *050

-126THROUGH LCOPF, FOR K1,l,K.G.M-1 *051 A(N+K) a -CF+CG/(1.+K*CH) *052 01 B(N+K) = 1./DELT+2.eCF *053 01 LOOPF C(N+K) - -CF-CG/(1.+K.CH) *054 01 A(N+M) a -CI *U55 B(N+M) - CJ+CI *056 ENTRYB THROUGH LCOPFA, FOR ItO,1,I.G.N *057 LOOPFA U(I) a INLET *058 01 THROUGH LCOPP, FOR J1,1,JJ.G.P *059 THROUGH LCOPG, FOR 1=0,1,I.G.N-1 *060 01 D(l) = F(I,J)/DELT+VEL(I)~CB*U(I)+Q 0061 02 LOOPG CONTINUE *062 02 D(N) * CODF(NJ)+0.5*Q/N *063 01 THROUGH LCOPH, FOR K=1,l,K.G.M-1 *064 01 LOOPH D(N+K) a W(KJ)/DELT *065 02 DIN+M) a CJ~W(MJ) *066 01 EPS(C) 8(0O) *067 01 GAMMA(O) = O(0)/EPS(O) *068 01 THROUGH LCOPK, FOR I~1,1,I.G.M+N *069 01 EPS(I) - B(I)-A(I)*C(I-1)/EPS(l-1) *070 02 LOOPK GAMMA(I) = (D(I)-A(I)*GAMMA(I-1))/EPS(I) *071 02 S(M+N) * GAMMA(M+N) *072 01 THROUGH LCOPL, FOR Is M+N-l,-1,I.L.0 *073 01 LOOPL S(l) z GAMMA(I)-C(I)*S(I+1)/EPS(I) *074 02 THROUGH LCOPM, FOR K-0,1,K.G.M *075 01 LOOPM VIK) - S(N+K) *076 02 THROUGH LOOPN, FOR I~0,1,I.G.N *077 01 U(I) a F(IJ) *078 02 LOOPN,(I,J)q S(I) *079 02 THROUGH LCOPO, FOR K-1,1,K.G.M *080 01 LOOPO W(K,J)~ V(K) *081 02 WHENEVER (COUNT+1)/FREQ.E.(COUNT+. )/FREQ *082 01 PRINT COMMENTSee.e**.*. eee-ee~eeweeee*.e*........~**.*. eSe *083 01 01 X * J/CB *084 01 01 PRINT RESULTS T,X *085 01 01 PRINT RESULTS S(O)...S(N)V(O)....V(M) *086 01 01 FGRAO- N-(S(N)-S(N-1)) *087 01 01 PRINT RESULTS FGRAD *088 01 01 WGRAD= M-(V(1)-V(O))/DELTA *089 01 01 PRINT RESULTS WGRAD *090 01 D0 CMPAREs WGRAD*LAM/BETA *091 01 01 PRINT RESULTS CMPARE *092 01 01 FACs 12.(1I.-RATIO).P.2./(RATIO.P.4.-4.*RATIO+3.)'093 01 01 SUMSt FAC*S(0)/8./CA *094 01 01 THROUGH LCOPS, FOR I-1,1,I.G.RATIO*N *095 01 01 LOOPS SUMS a SUMS+S(I)*FAC*I/CA *096 01 02 THROUGH LOOPT, FOR IIl1,I.G.N-1 *097 01 01 FAC- 12./(RATIO.P.4.-4.*RATIO+3.)-((1.-2.*RATIO)~I/CA+ *098 01 02 1 2.*RATIOII/(N.P.)-IP.3../N(N.P4.)) *098 LOOPT SUMS. SUMS+FAC*S(I) *099 01 02 MIXEDT- SUMS *100 01 01 PRINT COMMENT SOTHE MIXED MEAN TEMPERATURE ISS *1l -O' PRINT RESULTS MIXEDT *102 01 01 END OF CONDITIONAL *103 01 01 LOOPP CONTINUE'104 01 COUNT * CCUNT+1 *105 T = T+DELT *106 WHENEVER COUNT/FREQ.E.(COUNT+0.0)/FREQ *107 WHENEVER.ABS.((V(1)-V(D))*M/DELTA).G.ALPHA.OR.T.L.5.*X/U(0) *108 01 TRANSFER TO ENTRYB *109 02 END OF CONDITIONAL *110 02 TRANSFER TO START 11 01 END OF CONDITIONAL *112 01 TRANSFER TO ENTRYB i13 END OF PROGRAM 114

-127START. _____PREAD DATA - 1 PRINT DATA ^~r~~~\ ~ PRINT DATA f IS DATA Is DATA4 ~ REJECT STABLE F STATEMENT COMPUTE CONSTANTS INITIALIZE VARIABLES ITERATE DIFFERENCE EQUATIONS HAVE VARIABLES F CONVERGED - INCREMENT PRINT DATA _ RT REJECT T 1 IS CONVERGENCE STATEME NT ----- TOO SLOW F COMPUTE J>p PHIF AND AR V T COMPUTE | ^ PRINT WALL GRADIENTS RESULTS Figure 43. Computer Program B Flow Diagram, Steady-Periodic.

COMPUTER PROGRAM:! NOMENCLATURE STEADY-PERIODIC A(I) Sin rG coefficient in tpf at the point (I, J) in the fluid; A(I, J) ALPHA Convergence limit AMPRAF(I) Periodic amplitude ratio at the point (I, J) in the fluid, IA(I)2 + B(I)2 /e S(I) AMPRAW(K) Periodic amplitude ratio at the point (K, J) in the wall, c(K)2 + D(K)2 /e V(K) ATN1. Arctangent subroutine B(I) Cos rG coefficient in Vpf at the point (I, J) in the fluid; B(I, J) BETA Pf Cf/Pw Cw C(K) Sin rF coefficient in tpw at the point (K, J) in the wall; C(K, J) COUNT Iteration counter D(K) Cos rF coefficient in Vpw at the point (K, J) in the wall; D(K, J) -DELTA Dimensionless wall thickness, 6 E(I) A(I, J-l) EPS FGRAD Periodic fluid temperature gradient at wall FREQ Axial frequency of solution printout F(I) B(I, J-l) G(I) A(I, J) after the previous iteration. Used in the convergence inequality for A(I, J) GAM P i(3

-.129H(I) B(I, J) after the previous iteration. Used in the convergence inequality for B(I, J) I Radial counter in the fluid, (I = 0, 1, 2...N) J Axial counter, (J 0, 1, 2...P) K Radial counter in the wall, (K = 0, 1, 2...M) L Length of grid network in the axial direction LAM A M Number of radial elements in the wall MIXEDT Mixed mean temperature of fluid N Number of radial elements in the fluid P Number of axial elements in the wall and fluid PHIF(I) Of = Arc tan - B(I)/A(I) PHIW(K) =w - Arc tan - D(K)/C(K) PERCOF(I) Periodic coefficient in the fluid, A()2 +B(I)2 PERCOW(K) Periodic coefficient in the wall, C(K2 + D(K) Q(K) C(K, J) after the previous iteration RATIO Bingham plastic constant, a R(K) D(K, J) after the previous iteration S(I) tsf at the point (I, J), in the fluid, S(I, J) SQRT, Square root subroutine TIMEAV Time Average of periodic fluid temperature gradient, FGRAD U(I) S(I, J-l) VEL(I) Fluid velocity, VEL(:I) a u(r) VMAX Maximum fluid velocity, VivMX = uax/2 u* V(K) Vsw at the point (K, J) in the wall; V(K, J) W(I) S(I, J) after the previous iteration Y(K) V(K, J) after the previous iteration

-150COMPUTER PROGRAM B - STEADY-PERIODIC CIMENSION A(5),8B(5)),CI.(5C),O( 5L), F (5C) G ( 5C ) 0, 50 ),Q(50),R (5C), 5Cl 1 S > 5.J,L.1_,_E.[l5_ ) _ l5. )_LYI_-J.i, PHI W 1 5.), PH I F.(.5E ], PER COF ( 5 ), " 2 VEL(5),EPS(S(5),OPPF(50),CPPW50),ACJF(50),ACJh(5C), 01 3 AMPRAFt5T) wAPRAWO(50)PERCCW(5) 01 INTEGER I t J, K,F,N, P,COUNT,FREQ 02 PRINT COMPENT $2 HEAT TRANSFER TC BINGHAM PLASTIC IN CIRCULAR *'.03 1 TUBES IN LAMINAR FLOW WITH INTERNAL HEAT GENERATICNS C03 TA RT.... _.BEAL.. CATA *4 PRINT COPPENT $1NEW DATA$*905 PRINT RESLLTS EPSBETAtLAM,GAM,DELTA.,LP,RATIONM, FREO *06 PRINT COURENT $ THE NUPBER OF RADIAL PCINTS' IN THE FLUID IS$.r.07 PRINT REISLLTS N *008 CELT* 10.OE-37 -009 ALPHA =.C01*CQ1 PI = 3.14159 *0il WHENEVER RATIO.E.I. 012 THROUGH LCOPW, FCR I=O,lI.E.N *-13 01 LOCPW yEL(I)=').5'14 01 01 VY AX =0.5'.15 01 TRANSFER..TO EN.TRY,=16 Cl END CF CCNDITICNAL17 C1 VMAX = 3.. 1.-RATIC).P.2.)/.(RATIC.P.4.-4.*RATIC+3.) *018 THROLGH LCCPX, FCR IC,t1I.IE.N *019 IF=I *020 01 KNFN *021 01 _IIELLL-. 3.-iL-2ll2. BA.TI CZ.RAT.L,.LELNF- IF.* LEIE LNFIy. *022 01 1 RATIC.P.4.-4.eRAT(C+3.) *022 tCCPx WHHENE.VER IE/NF.LERA.TIOi VELII)=.VMAX *.23 01 ENTRY THROUGH LCOPT, FOR P=P,5,.P/L.G.5C00~ *024 WHENEVER P/L.G.GAM/VELIN-1) *025 01 PRINT COPFENT BOTHE NUPBER OF AXIAL NCCAL POINTS IS.$ *026 01 01 ER INI_..RE SLLI&S_..P t~02.7..1 01 TRANSFER TO ENTRYA *028 01 01 CThERWISE *029 01 01 CONTINUE *!30 01 01 LCEp.T END CF CCNDITIONAL *031 01 01 PRINT COVPENT $CCIVERGENCE OCCURS IN THE FLUIC EQUbITIONS WHEN *032 1_P.PIS INCRFASIDn Tn P/LI. 5DnL CATA REJECTE.032 TRANSFER TO START *C33 EN-RYA IHROUIGH.LCOPL, FCR..M*M -M.L2.....*0.34 WHENEVER GA.G. (2.*N-l.+2.*LAM*//BETA/CELTA+LIM/BETA) / OELTA/ *035 01 1.. BETIA12. _. t,0 5./N.1.AND,.GAM G..4. *LA _.p,*M./~ ELAY D.ELT A..035.. PRINJ COPPENT $0 THE NUMBER OF NCDAL POINTS IN THE WALL,hb *036 01 01 PRINT REsI;.LTf. *037 01. 01 TRANSFER TO'ENTRYB *038 01 01 CTHERWISE *Q39 01 01 CONTINUE *040 01 01 LECPtM -.EJA. CE.CCN.IOINAL- *41 01 01 PRINT COPPENT$O CIVERGENCE OCCURS IN THE WALL EQUATION6 WHEN *042 1 M IS RFPtEC.F To 2. DATA RFJFCTEDI **042 TRANSF'ER TO START *043 C8 = 1.*P/L *045

-131 CC = CO/GAM "'46 CC = I/CELTA+3.5 *C47 CE_ LAMIGAY *48 CF P= FM/CEl TA/CELTA *2149 CG GAM*CELTA/LAM/2./~ *.54 CH LAMe(M/CELTA+M-;.5) -051 CI = GAMeCEITA*3.5/ Me(1.+CELTA) *. 52 PRINT COPPENT $:THE CONSTAKTS GENERATEC eY T-IS CATA ARES *~53 PR.IN.T_ RE.SLLTS CACE,C,CC,CC,CE,CFCG,C -,C I,VMAX, N EL ( C )..VEL I N-1) *-54 THRCtGH LCOFPA, FCR I[C,~1I.G.N *55 (I) =. C',56 C1 B(I) =. *.57 01 S(I) = C. ~*5R 01 EUL) = 2. *059 01 LGGEPA LI) * C.061 C1 THROLGH LCCFV, FCR I'l0,II.G.M *062 C( ) =. *('63 01 LCOPV CII). *064 01 THROLGH LCCPE, FCR Jtlll,J.G.P *665 T.HROLGH LCCPC. FOR I0~l,,I.G.N *066 01 G(I = 0. *?67 02 I- ( = 0. *068 C2 LCOPC h Il) = 0. *C69 02 THRCLGH LCOPC, FOR I*0,1,I.G.M *070 01 CU() = 0. *971 02 LC0OCE R_.LL.__... -072 C2 CCUNT= 0 ~073 01.EhRYC C OM)' CH*(CIP)-C(M —1))/CI *074 01 C(M)= -CI-eLU(M)-C(M-l) )/CI *075 01..THROUGH LCOPE, FCOR K"M-1t,-lK.L.1 *076 01 FAC= 2.*CELTA*.(KODELTA/M+. ) *077 02 C(K ) -CE* CFf<M/FAC)*C (K+l4]-CE*(CF-Y/FAC C(K-l) *078 02 1 +2.*CECF*Cr.(K) *078 ~LCPEe.'E C.(.K)*. C,.(CF+l/FAC)*0(K+I)+CE*.(CF-M/FAC)*DtK-1)-2.4CE*CF*0(K) *079 02 BIN)= L-LAM/BETA~CD'.IC.(l)-C,(IO)+(N-0.5).(AN)-A(N-l) )-EPS *080 01 1 12,/K) / )tGAMH (OELTA/BETA/2/M0+5IN ] 2 ~80 A(N)~= LAM/HETA*CDO(DI l)-D(O) )-(N-C.5)(B(N)-BIN-1) ))/ ~081 01 1 GAMI ( CELTA/LEETA/2./M+0.5/N) *'081 SIN)' 0.5/ytNN-0.5N) +S(N-1) *082 01 C.O)' AMI) *083 01 CIO) B(N) *084 01.THROUGH LCOPFi FOR /'N-l t1tl-.LL *085 01 Ar~ ) (CA*ll.+C.5/I )cA(+1 )*CA( 1.-0.5/I )A.Jl -l)*CBe.VEL Cl) *086 02 1 EL ),+GAMe*!I) +EPS)/ 2.eCA+CB#VEL ( I) 086 e.{)1 )CA I1.+0.5/I)*B(I+1)+CA*( l.-0.5/I)*B(I-l)+CB~,VEt(:() *087 02._ FI I )-.GA/. A(.I ).1L (Z....C A+..*y_EL...).) 087. LBPF S$l)~= kCA 1.*0.5/,I)*S(I+1 )CA.(l.-C.5/1)S.(I-,l) +C8>.U(I) *088 02'1....vE.(. 1+L.jt1 LZ.)_! CA +CB6VE L iI))) 1 -088 A10)k t4.ECA'A(1l)+C8,VMAX*E(O)+GAMRBIO)+EPS)/ *089 01 1 4, *CA+CB.VMAX *089 8(0)' 14.'*CA'B.(1)+CBVMAX*FIU)-GAMAlIUJ 1)/ *090 01'.[lCCi+CBeV.8'AX}) *090 St(O) t44.CA,'S( 1)+CB6VMAX*U(0O)+1l. )/(CBeVMAX+4.eCA)'091 01 THRO GH. LC P G..FR.. O, L,..G, M. 0_92 01 NHENEVER.ARS...(C(l))L.DELT.OR-.AB6.(0lI )iL.CDELT *093 02 TRANSFER To ENTRYE~ 094 01 02 CTHERWISE *095 01 02 IHEN.VER..AP lS.,ICiOZJI: IIJ_.G.,I0OLCOflANll..._ABF&,CLLYOIlL.-fI lCGOfOC *096 01 02 TRANSFER TC LCOPG *097 02 02

-132ENC CF CCNDITICNAL *r.98 22'2 WHENEVER.APS. ((C( I )-CG(I) )/C( I)).L.ALPf-I9..ANC.:99 1 C2 1.ABS-( LD(Il-P ( IJ)/C( I)).L.ALPHA~1 *~C99 CCNTINUE *l'(* 02 C 2 CTHERWISE.1* 02 r2 CCUNT=CCLNT+1 -132 02 C2 WHENEVER CCUNT.G.lC' -13 f,2 2 PRINT CCFFENT 4***4********14*~*< <* * *U*4 E *4I* ****.$ *1~4 03 02 PRINT COPFENTs CCNVERGENCE IS TCC SLCW. DATA REJECTEC$ 105 3 02 PRINT CCPENKT$ SLOW CCNVERGENCE IN THE WALL ECLATICNS.$'*16 3 02 PRINT RESLLTS S ( )...S(N),W( )...W N),A(1)...A), G (),...G(N), *1r7'3 C2 1 f8(0)).. B(N).K (O)...H(N),C(<;)...*C()~)~C(f).. oC(M) t~107 2 R(0)...R(M) ~ -1 7 PRINT COPfENT $**~* * * ~*******U*~~*~ 4t*~*~~~*~~~~~ $* ~1,R 03'2 TRANSFER TC. START -139 3.12 CTHERWISE *110 C3 02 CONTINUE *11l1 3 02 END CF CCODITICNAL *112 C3 02 THROLGH LCCFH, FCR K=1,1.K.G.M *113 02 12 C(K) = C(K) *114 02 03 LCCEI- R(K = D(K) *115 02 33 THROUGh LCCPI, FCR K*C~l,K.G.N *116 02 2 GIK) = A(K) *117 C2 C 3 H(K) = BE(K) *118 C2 23 LCEF.I WtK) = S(K).*119 02 03 TRANSFER TC ENTRYC.121 02 T2 END CF CCNCI.TICNAL *121 C2 02 END CF CCNDITICNAL *122 01 02 LCCPG CONTINUE *123 02 THROUGH LCCPJ~ FOR If=*l,I.G.N *124 l1 WHENEVER.AdS.(A(I)).L.DELT.CR..ABS.(B(I)).L.CELT.CR..ABS. *125 02 1 (S(I)).L.CELT *125 IRALNSFER JO ENTRYE.. *1........ -26 01 02 OTHERWISE *127 01 q2 WHENEVER.AHS.I(A(I)-G(I))/A(I)).L.ALPHA.ANC;.A8S.((B(I)-H(I))/ *128 01 02 1 Etl)).L.ALPA.ANC....ABS.(SlI)-W(I))/S(I)).L.ALPHA *128 CONTINUE *129 02 02 CTHERWISE *130 02 C2 CC U.N ICOL;T+1. _.._. *131 02 02 WHENEVER COhNT.G.25.ANC..ABS.(B(I)).G..ABS.(1.0E C5A(I)) *132 02 02 PRINT COPENIT $*133 03 02 PRINT COPFENT $SLOW CONVERGENCE IN FLUIC EQUATIONS.EECALSE A *134 03 32 1 TENCS TCWARC ZEROS *134 PRINT RESLLTS W(C)...W(NfGfGf...G(N),H(O)...H(N) *135 03 02.IRANSEER10-_ENLRY. *136 03 02 END OF CCKDITICNAL *137 03 02 WHENEVER CCLKT.G.lCO *138 02 02 PRINT COMFENT $ ***o*ooeeo~e~ eeremeee, ni.e 0 emoee eee$ *~139 03 2 PRINT CORNEN-f$ CCNVERGENCE IS TOO SLOW. DATA REJECTEDS *140 03 02 PRINT COFENTS SLOW CONVERGENCE IN FLUIC EQUATIONS.$ *141 03 02..e. LN[__F S L T S S (I *_.CIL. LW.I), A (o).l.iA. L...G.1... I, J *142 3 02 8(f)...B(N).H(O)...H(N)C(0~)...C(M),C(0)...C(),C(C)...D(M), "~142 2 RIO)...-R_(..1 *142 PRINT COFPENT $~<i~m~ffe* ~~, **mmee*,,,m~4 ~ *04~~ 0~~e~m~0e~~$ *143 03 02.TRANSFER _J._ART.*144 03 02 OTHERWISE *145 03 02 rnNTZINjEL 4603 02 END OF CCKCITICNAL *147 03 02 -THkNg n H-tG OIR FaLY4_FR K*A"_l*fK.GMJ *148 02 02 CIK) f C'(K) *149 02 03

-133LCCPK RIK) = D(K) *15 202 C3 THROUGH LCOPL, FCR K0t,1,K.G.N *151:2 C2....G[.._..... _______ ~*K)- A iK.)_.152 02 "?3 H(K) = B.(K) *153 02 53 LO3QPL WAiK) = SI(K) 154 02 C3 TRANSFER TO ENTRYC -155 (2 "2 END CF CONDITIONAL'156 C2 ('2 EhTRYE END CF CCNCITICNAL *157 01 02.... -LELd; ONrT_[.-;.l.........._.._........ * 159 C2 EITRYD THROLGH LCCP?', FCR I=C,1I.G.N *159 01 Elt). A.(I) 169 2 F{{) = B(I)'161 02 *LC.CfiP.... L(I) = S(1)'162'2 WHENEVER J/FPRE.E.(0.G+J)/(0O.OFREC)'163 C1 _-...............x. ~=JCB._..._....__.___............................'*164 01 01 PRINT COO!ENT $0$ *165 01 Cl PRINT RESLLTS X,CCUNT *166 01,1 THROUGH LCChN', FOR IO',1,I.G.NM 167 C1 71 NCRMF 2 l./(.ABS.(8(I))) *168 Cl 12 CFPF(I) =,-R(I)*NORMF *169 01 02 ADJFL._A(I)*NORMF _ 170 01 32 PHIF(I) = ATNK.(CPPF(I)oADJF(I)) *171 01,02 tCCPIK PERCOF(I) = (SQRT.((CPPF(I).P.2.)+(ACJF(I).P.2.)))/KCRMF *172 O1 02 THROUGH LCOPC4 FCR KfO,l,lK.G.M *173 01 01. NORMW a l./I.ABS.(D(K))) *174 01 02 CFPWIK) = -~(K)*NORMW *175 11 l C2 _ ACJW =CLK__ )NORMW ___.....176 f 1:02 PHIWIK) = ATI. (CPPW(K) ADJW(K)) *177 01 02 LC...... PERCOW(K)'" (SCRT.((OPPW'(K).P.2.)+IACJWIK).P.2.)))/NCRMW *178 01 C2 THROtGH LCCPC, FCR IOt,1,I.G.N *179 C1 1 L....... AMPRAF(I) " PERCCF(I) / (S ( I)*EPS) -180 Cl1 02 THROUGH LCOPR-, FOR I=0,1,I.G.M *181 C1 31 LCEPR AMPR AW(I) U.= CCW ( I/S (N_)_EPS). *182 01 2 FAC= 12.(lI.-RATIO).P.2./(RATIO.P.4.-4.*RATIC+3.) *183 01 01... SUMS' FAC3S~C)/8./CA'184 01 01 SUMA= FAC*AIC)/8./CA *185 01 11..SliMe FAC*I8C)/8./CA *186 01 C1 THROIGH LCO;Pi FOR It1l,l,I.G.RATIOGN *187 01 01 SUMS = SUS+StI)JFACI/C_ A.. -188 01 02 SUMA = SUPA+/AI)eFAC*I/CA *189 01 C2..J...IOP SUMB = SLFB+e(II)eFACCI/CA'190 01 C2 THROUGH LCOPZ, FOR I1I,1,I.G.N-1 *191 Cl C1 FAC= 12. / RATIO.P.4.-4.*RATIO+3.) ( 1.-2.RATIC )*I/CA+ *192 01 02 1 2.*RATIOI* I/(' N.P.3. )I-I.P.3./(N.P.4.)) *192 SUMS= SUMS~+_FAC*S( I )..____._....__...___... - 193 l1 02 StUM:A SUMA+F0ACA(II) 194 01 02 ICP.Z............SUMEB' SUM8*FAC8( I ) 195 01 02 NUMA= Ne(AIN) —A(N-1)) *196 01 i1... -NMB= N(BILN)-B(N-l))'197 (1 01 PRINT COIPENT SOTHE STEADY PERIOCIC WALL TEMPERATURE *198 01 01 i GRACIENT IN PI/6 INCREMENTS ISL __$ _._....................... -_ 198 SUM= C.'199 01 01.-... T...H..ROLGH LCCORS, FOR ARG=PI/6.,,PI/6.,.ARG.G.2'*PI *2CC 01 01 FGRAC= NU)A*SIAN. U (ARG) N201 01 02... SUM= SUM+FGRAD'202 01 02 MIXECT='SUMS+SUMA*SIN.(ARG),+SUMBCOGS.(ARG) 203 01 02 lCEPS_ PRINT BESL LTS ARG mFGRAC,'M XE1O.. — T...... -20)4 01 02 TIM'EAV* SIM/12. 205 01 01... RINJ.COJ''EN.T..OTHE.TIME AVERAGE OF THE WALL TEMPERATURE'2n6 01 01 I GRADIENT ISS..206....... _P.E.T.REfSL. Sl T 5.TiMEAV......................... 2C7 1 01 PRINT RESILTS S(0)...S(N)hPHIFIO)..PHIF(N);tAMERAFIO)...AMPRAF(N) 208 01 01 1 AaIO A IN I.... lIN).PHI WIOQ)...eL )............... 208 2 AMPRAW1[)...AMPRAW(M).CLO),..CIM),D(0)...CLM) *208 _D....... __.___.(4D F CCNDII{ONAL....2C9 01 01 LUCCO1 CONTINUE *21,? 01 _..IR NSFER. TO. ITART -211 END CF PROGRAM'212

BIBLIOGRAPHY 1. Schechter, R. S., and E. H. Wissler, "Heat Transfer to Bingham Plastics in Laminar Flow through Circular Tubes with Internal Heat Generation," Nucl. Sci. and Eng. 6 (1959) 371-375. 2. Poppendiek, H. F., "Forced-Convection Heat Transfer in Pipes with Volume-Heat Sources within the Fluids," Chem. Eng. Sympos. Ser., 50 (1954) 93-104. 35 Grigull, U. "Warmeubergang an nicht-Newtonsche Flussigkeiten bei laminarer Rohrstromung," Chemie-Ing. Techn. 28 (1956) 553-556. 4. Sparrow, E. M., and R. Siegel, "Laminar Tube Flow with Arbitrary Internal Heat Sources and Wall Heat Transfer," Nucl. Sci. and Eng., 4 (1958) 239-245. 5. Tachibana, F., and T. Morishita, Heat Transfer of Inside Tube Flow of Slurry, Industrial Research, Tokyo Univ., 11 (1959) 45-46. 6. Michiyoshi, I., Heat Transfer of Slurry Flow with Internal Heat Generation - Part I. Bulletin of Japan Soc. Mech Mech. Eng. 5 (1962) 315-319. 7. Michiyoshi, I., et al., "Heat Transfer of Slurry Flow with Internal Heat Generation - Part II," Trans., Japan Soc. Mech. Eng., 28 (1962). 8. Ayers, D. L.,'The Asymptotic Behavior of Periodic Heat Transfer to Laminar Flow in Tubes, Ph.D. thesis, Department of Mechanical Engineering, University of Michigan, (1963). 9. Clark, J. A., V. S. Arpaci and K. M. Treadwell, "Dynamic Response of Heat Exchangers Having Internal Heat Sources - Part I", Trans. ASME, 80 (1958) 612-624. 10. Arpaci, V. S., and J. A. Clark, "Dynamic Response of Heat Exchangers Having Internal Heat Sources - Part II." Trans. ASME, J. 80 (1958) 625-634. 11. Arpaci, V. S., and J. A. Clark, "Dynamic Response of Heat Exchangers Having Internal Heat Sources - Part III," Trans. ASME, J. of Heat Transfer, 81 (1959), 253-266. 12. Yang, W. J., J. A. Clark, and V, S. Arpaci, "Dynamic Response of Heat Exchangers Having Internal Heat Sources - Part IV," Trans. ASME, J. of Heat Transfer, 83 (1961), 321-323. -154

-13513. Yang, W. J., "Dynamic Response and Resonance Phenomenon of a Single-Solid, Single-Fluid Heat Exchanger - Part I," Trans. Japan Soc. Mech. Engrs. 27 (1961), 1276-1285. 14. Yang, W. J. "Dynamic Response and Resonance Phenomenon of a Single-Solid, Single-Fluid Heat Exchanger - Part II," Trans. Japan Soc. Mech. Engrs., 27 (1961), 1892-1907. 15. Yang, W. J. "Dynamic Response and Resonance Phenomenon of a SingleSolid, Single-Fluid Heat Exchanger - Part III," Trans. Japan Soc. Mech. Engrs. 28 (1962), 551-558. 16. Keller, R. B., The Dynamic Response of a Heat Exchanger with Constant Internal Heat Generation to Several Types of Fluid Velocity Transients," Ph.D. thesis, Department of Mechanical Engineering, University of Michigan, (1961). 17. Yang, W. J., "Frequency Response of Heat Exchangers Having Sinusoidal Internal Heat Generation," ASME paper 62-H-21 presented at the 5th National Heat Transfer Conference, ASME-AIChE, Houston, Texas, (1962). 18. Yang, W. J. "Temperature Response of Nuclear Reactors Having Sinusoidal Space-Dependent Internal Heat Generation," 1961 National Nuclear Congress, Abstract No. 48, Tokyo, Japan and Trans. Japan Soc. Mech. Engrs., 28 (1962). 19. Metzner, A. B., "Non-Newtonian Technology: Fluid Mechanics, Mixing, and Heat Transfer," Advances in Chemical Engineering, Volume 1, Academic Press Inc., New York (1956). 20. Bird, R. B., et al., Transport Phenomena. John Wiley and Sons, Inc., New York (190j). 21. Yang, W. J., Transient Heat Transfer in Heat Exchangers Having Arbitrary Time- and Space-Dependent Internal Heat Generation, Publication of the Industry Program, College of Engineering, the University of Michigan, IP 582 (1962). 22. Harnett, J. P., "Experimental Determination of the Thermal Entrance Length for the Flow of Water and of Oil in Circular Pipes," Trans. ASME, 77 (1955), 1211-1220. 23. Eckert, R. G., and A. J. Diaguila, "Convective Heat Transfer for Mixed, Free, and Forced Flow through Tubes," Trans. ASME, 76 (1954) 497-504.

-13621. Schneider, P. J., "Effect of Axial Fluid Conduction on Heat Transfer in the Entrance Regions of Parallel Plates and Tubes," Trans. ASME, 79 (1957), 765. 25. Lane, J. A., H. G. MacPherson, and Frank Maslon, Fluid Fuel Reactors, Addison-Wesley Publishing Co., Inc., Reading, Massachusetts, (1958). 26. Carnahan, B., H. A. Luther, and J. 0. Wilkes, Applied Numerical Methods, John Wiley and Sons, Inc., New York, (1964). 27. Metais, B., and E. R. G. Eckert, "Forced, Mixed and Free Convection Regimes," Trans. ASME, J. of Heat Transfer 86 (1964), 295-296. 28. Thomas, D. G., E. L. Compere, and J. P. McBride, "Thorium Oxide Suspensions," Nucleonics, 18 (1960), 104-110. 29. Thomas, D. G., "Homogeneous Reactor Project Quarterly Progress Report for the Period Ending July 31, 1954," USAEC Report ORNL-1813 (Del.), Oak Ridge National Laboratory, (1954), 125. 30. Beck, J. V., "Thermocouple Temperature Disturbance in Low Conductivity Materials," Trans. ASME, J. of Heat Transfer, 84 (1962), 124-132. 31. Inman, R. M., "Experimental Study of Temperature Distribution in Laminar Tube Flow of a Fluid with Internal Heat Generation," International Journal of Heat and Mass Transfer, 5 (1962), 1053-1058. 32. Pauling, L., College Chemistry, W. H. Freeman and Co., San Francisco, (1956). 33. Pannell, J. R., "Experiments with a Tilting Manometer for Measurement of Small Pressure Differences," Engineering (Sept. 12, 1963) 343-344. 34. Tachibana, F., T. Morishita, T. Naito and M. Okubo, "Pool Boiling Heat Transfer in Slurries," Paper no. 28 presented at the Heat Transfer and Thermodynamics Conference, Japan Society of Mechanical Engineers, (Nov. 28, 1963), Nagasaki, Japan. 35. Orr, C. Jr., and J. M. Dalla Valle, "Heat-Transfer Properties of Liquid-Solid Suspensions," Chem. Eng. Progr. Symposium Series no. 9, (1954), 29. 36. Streeter, V. L., Fluid Mechanics, Third edition, McGraw-Hill Book Co., Inc., New York (19627). 37. Kline, S. J., and F. A. McClintock, "Describing Uncertainties in Single-Sample Experiments," Mechanical Engineering, 75 (1953), 3-8.