THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE ASYMPTOTIC BEHAVIOR OF PERIODIC HEAT TRANSFER TO LAMINAR FLOW IN TUBES David L. Ayers 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 1963 May, 1963 IP-622

I

PREFACE The author gratefully acknowledges all those who aided in this work and in particular the following: Professor V. S. Arpaci, Chairman of the author's Doctoral Committee, for his invaluable advice and encouragement throughout the course of this worko Professor J. A. Clark, A. G. Hansen, S. W. Churchill, and Go J. Van Wylen, members of the Doctoral Committee, for their interest and cooperation. Mr. Brice Carnahan, Assistant Director, Ford FoundationProject on Computers, for advice on computer programming. The Department of Mechanical Engineering for providing financial aid in support of the experimental investigation. The Computing Center, Dr. R. C. F. Bartels, Director, for providing gratuitous computer time. The Horace H. Rackham School of Graduate Studies for financial aid in the form of fellowships. The Industry Program of the College of Engineering for aid in the preparation of the manuscript. 111

TABLE OF CONTENTS Page_ PREFACE..oooo. ooo....................................... iii LIST OF TABLES....,...,.................o o... vi LIST OF FIGURESo a..... oa.. oo..................a.a........... vii NOMENCLATURE.o.o.o....o....o......a..*........*.........*..... x COMPUTER PROGRAM NOMENCLATUREo.o........... o o o............ xv CHAPTER I INTRODUCTION......................... 1 II ANALYSIS........OO...............o.eo.o. ee. 4 A. Introduction...................................... 4 B. Formulation8...... 8 C. Finite Difference Equations....................... 12.D. Computer Solutions..., 14 1. Stability Tests................... 18 2. Convergence Tests........................ 26 3. Nusselt Number Computation.................... 28 III NUMERICAL RESULTS................................. 34 A. Steady Solution................................ 35 B. Steady-Periodic Solution....................... 45 IV EXPERIMENTS............ oooo................. 68 Ao Apparatus................... 68 1. Test Section................................. o 68 2. Fluid Supply System. o.... o..*..,0......... 72 3. Line Heater.........................o........ 72 4. Power Supply........ o........................ 73 B. Instrumentation.......O...... *.................... 73 1. Thermocouples................................ 77 2. Voltage and Current........................... 80 3. Temperature Recording.......................... 80 4. Mass Flow Rate............................. 82 iv

TABLE OF CONTENTS (CONT'D) Page C. Test Procedure...........-.................. 82 D. Data Reduction...................................... 83 Eo Results of Experiments............................. 86 V CONCLUSIONS.......................o........ 91 APPENDICES..........................,................... 93 A DERIVATION OF THE DIFFERENTIAL EQUATIONS................ 93 B DERIVATION OF THE FINITE DIFFERENCE EQUATIONS........... 97 C COMPUTER PROGRAMS.................................... 118 D THIN WALL TEMPERATURE DROP.............................. 127 E THE EFFECT OF AXIAL CONDUCTION......................... 134 BIBLIOGRAPHY............................................. 136 V

LIST OF TABLES Table Page I Difference Equations - Large r..................... 15 II Stability Inequalities - Large r Equations.......... 22 III Difference Equations - Small r.................... 2 vi

LIST OF FIGURES Figure Page 1 Schematic Diagram of Physical Model................... 2 Free Convection in a Fluid in a Horizontal Heated Tube...............................e. 7 3 Effect of Free Convection on a Normal Laminar Temperature Profile in the Entrance Region............. 7 4 Concentric Grid Network......................... 16 5 Steady Nusselt Number versus Axial Distance............ 36 6 Periodic Nusselt Number versus Time........37 7 Steady Temperature at r = 1.0 versus Axial Distance for Various Wall Thicknesses..................... 39 8 Steady Fluid Temperature versus Radius for Various Axial Distances;;The Development of the Steady Temperature Profile in the Fluid..................... 40 9 Steady Wall Temperature versus Radius for Various Axial Distances; The Development of the Steady Temperature Profile in the Wall. 42 10 Steady Temperature at r = 1.0 versus Axial Distance for Various /A............................. 43 11 Steady Wall Temperature versus Radius for Various P/A................................44 12 Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r = 1.0.47 13 Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r = 0.9........ 48 14 Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r = 0.8........ 49 15 Fluid Temperature Amplitude Ratio versus Axial Distance for Various Wall Thicknesses, r = 1.0.... 51 16 Fluid Temperature Phase Angle versus Axial Distance for Various Generation Frequencies............ 53 vii

LIST OF FIGURES (CONT'D) Figure Page 17 Fluid Temperature Phase Angle versus Axial Distance.... 55 18 Fluid Temperature Amplitude Ratio versus Axial Distance for Various Generation Frequencies............ 56 19 Fluid Temperature Phase Angle versus Axial Distance for Various Radii, r = 100................... 57 20 Fluid Temperature Phase Angle versus Axial Distance for Various Radii, r = 500.........0...... 59 21 Fluid Temperature Phase Angle versus Axial Distance for Various Radii, r = 1500................. 60 22 Fluid Temperature Amplitude Ratio versus Axial Distance for Various Radii, r = 100 e0*-0o..- 61 23 Fluid Temperature Phase Angle versus Radius for Various Axial Distances......,..,.................. 62 24 Fluid Temperature Amplitude Ratio versus Radius for Various Axial Distances............................ 64 25 Fluid Temperature Phase Angle versus Generation Frequency for Various Axial Distances.................. 66 26 Fluid Temperature Amplitude Ratio versus Generation Frequency for Various Axial Distances........... 67 27 Experimental Apparatus Showing Test Section........... 69 28 Experimental Apparatus Showing Reservoir Tower and Controls...........-...... on0o. 70 29 Schematic Diagram of Experimental Apparatus........... 71 30 Heater Exit Fluid Mixing Baffle........... 74 31 Experimental Apparatus Showing Four-Bar Linkage and Power Supply Auto-Transformer.............o..... 75 32 Schematic Diagram of Instrumentation...... 76 33 Typical Installation of Wall Thermocouples............ 78 viii

LIST OF FIGURES (CONT'D) Figure Page 34 Thermocouple DC Pickup Calibration............... 79 35 Typical Visicorder Trace during SteadyPeriodic Operation................................... 84 36 Steady Wall Temperature at r = Ro/Ri versus Axial Distance........................................ 87 37 Schematic Diagram of the Actual Wall Temperature Distribution at the Entrance of the Tube.............. 88 38 Steady-Periodic Wall Temperature at r = Ro/Ri versus Time...................................... 90 39 Control Volume in the Fluid........................... 94 40 System in the Wall.................................... 94 41 Wall and Fluid Grid Network......................... 98 42 Fluid-Wall Interface Grid Element.........o.......... 109 43 Grid Element at the Outside of the Wall............... 112 44 Fluid-Wall Interface Grid Element for a Radially Lumped Wall............................. 115 45 Computer Program Flow Diagram.......................... 119 ix

NOMENCLATURE A Coefficient of Sin ct in the periodic part of the dimensionless temperature in the fluid Ay 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; Aw, area of an element in the wall; Af, area of an element in the fluid; ft2 AR Amplitude ratio; ARf = JAs2+B2/esf; ARw =- C2+D2/ Esw; AR*, measured amplitude ratio at the outside of the tube; dimensionless B Coefficient of Cos ct in the periodic part of the dimensionless temperature in the fluid C Coefficient of Sin ct in the periodic part of the dimensionless temperature in the wall Cf Fluid specific heat, BTU/lbm-OF Cw Wall specific heat, BTU/lbm-OF D Coefficient of Cos ct in the periodic part of the dimensionless temperature in the wall Di Inside diameter of tube, ft g Local acceleration of gravity, ft/sec2 h Coefficient of heat transfer between the wall and fluid; hp, periodic heat transfer coefficient defined in Equation (II-19); hp, time average of the periodic heat transfer coefficient; BTU/hr-ft2_~F i I { } Imaginary part of the complex function in brackets kf Fluid thermal conductivity, BTU/hr-ft-OF kw Wall thermal conductivity, BTU/hr-ft-OF L Length of test section, ft.'m Fluid mass flow rate, lbm/hr

Nu Nusselt number, hDi/kf; Nus, steady Nusselt number defined in Equation (II-27); Nup, periodic Nusselt number defined in Equation (II-22); -Nup, time average of the periodic Nusselt number; dimensionless Pi Radially dependent coefficient defined in Equation (D-ll) P2 Radially dependent coefficient defined in Equation (D-ll) P2 Radially dependent coefficient defined in Equation (D-14) P4 Radially dependent coefficient defined in Equation (D-14) P5 Radially dependent coefficient defined in Equation (D-17) P6 Radially dependent coefficient defined in Equation (D-14) p Pressure, lbf/ft2 Pr Prandtl number,.ICf/kf, dimensionless q" Volumetric heat generation rate, BTU/ft3-hr qo Mean volumetric heat generation rate, BTU/ft3-hr R A specific radius; Rf, 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; ft R { } Real part of the complex function in brackets pUDi dmensionless Re - Reynolds number, dimensionless r Dimensionless radius, r*/Ri r* Radius, ft T Temperature; Tf, fluid temperature; To, entrance temperature; Tw, wall temperature; Tmm, mixed mean temperature, ~F t Time, sec U Fluid velocity; Um, average fluid velocity; Umax, maximum fluid velocity; ft/sec x Dimensionless axial distance x* Axial distance, ft y Dummy variable, dimensionless xi

Z Complex dimensionless temperature Q + itpw, defined in Equation (D-6) caf Fluid thermal diffusivity, ft2/hr aw Wall thermal diffusivity, ft2/hr Pfi Dimensionless parameter, pfCf/PwCw r Dimensionless frequency, R2i/af Dimensionless wall thickness, Ro 1 Ri bf Thickness of an element in the fluid, ft 5w Thickness of an element in the wall, ft E Dimensionless maximum periodic heat generation amplitude, defined in Equation (A-l) Dimensionless parameter, 5. r e Dimensionless complex temperature coefficient defined in Equation (D-7) Dimensionless time, taf/R2 A Dimensionless parameter, aw/af Fluid viscosity lbm/ft-sec of Fluid density, lbm/ft3' PW Wall density, lbm/ft3 a Dimensionless parameter, VIT. Ro A Ri Phase angle defined in Equation (D-14) Phase lag of temperature with respect to the heat generation function; of, phase lag of the periodic dimensionless temperature in the fluid; Ow, phase lag of the periodic dimensionless temperature in the wall; 0, measured temperature phase lag at r = Ro/Ri; radians cp Phase angle defined after Equation (D-17) 4t ~ Dimensionless temperature xii

4f Dimensionless fluid temperature'rWy Dimensionless wall temperature ~* Measured dimensionless wall temperature at r = Ro/Ri lpf Periodic part of *f 4pw Periodic part of *w *sf Steady part of *f *sw Steady part of *w At* Measured steady dimensionless temperature at the outside of the wall Measured periodic dimensionless temperature at the outside of the wall IQ Dimensionless wall temperature analogous to *pw but with a Cos wt generation rate Xu Frequency, radians/sec Bessel Functions Io[y] Modified Bessel function of the first kind, zeroth order, of the variable yo Il[Y] Modified Bessel function of the first kind, first order, of the variable y. Ko[y] Modified Bessel function of the second kind, zeroth order, of the variable yo K1[Y] Modified Bessel function of the second kind, first order, of the variable yo BER[y] R {Io[iJy]} BERl[y] - {Il[i2y} BEI[y] I {io[iJy]} BEI1[y] R fIl[ify]} KER[y] R {Ko[T~y]} xiii

KER1[Y] I {Kl[i y]} KEI[y] I {Ko[i ~y]} KEIi[y] -R {K1[i2 y]} (Reference 17 gives the details of the relationships which exist between the Bessel functions used here.) xiv

COMPUTER PROGRAM NOMENCLATURE A(I) Sin F9 coefficient in Vpf at the point (I,J) in the fluid; A(IJ) ALPHA Convergence limit defined in Equation (II-22) AMPRAF(I) Periodic amplitude ratio at the point (I,J) in the fluid, IA(I)2 + B(I)2/eS(I) AMPRAW(K) Periodic amplitude ratio at the point (K,J) in the wall,'C(K)2 + D(K)2/eV(K) ATN1.o Arctangent subroutine B(I) Cos FP coefficient in 4pf at the point (I,J) in the fluid; B(I,J) BETA PfCf/PwCw C(K) Sin rG coefficient in 4pw at the point (K,J) in the wall; C(KJ) COUNT Iteration counter D(K) Cos PF coefficient in Vpw at the point (K,J) in the wall; D(K,J) DELTA Dimensionless wall thickness 5 E(I) A(I, J-1) EPSQ C F(I) B(I, J-1) G(I) A(I,J) after the previous iteration. Used in the convergence inequality for A(IJ) GAM r H(I) B(I,J) after the previous iteration. Used in the convergence inequality for B(IJ) I Radial counter in the fluid, (I = 0,1,2...N) J Axial counter, (J = 0,1,2.o.P) K Radial counter in the wall, (K = 0,1,2.,.M) XV

L Length of grid network in the axial direction LAM A M Number of radial elements in the wall N Number of radial elements in the fluid NUSTED Steady Nusselt number, Nus NUSPER Periodic Nusselt number, Nup P Number of axial elements in the wall and fluid PHIF(I) Of = tan-1 -B(I)/A(I) PHIW(K) Ow = tan-1 -D(K)/C(K) PERCOF(I) Periodic coefficient in the fluid, IA(I)2+B(I)2 PERCOW(K) Periodic coefficient in the wall, VC(K)2+D(K)2 Q(K) C(K,J) after the previous iteration R(K) D(K,J) after the previous iteration S(I) *sf at the point (I,J), in the fluid, S(I,J) SQRT. Square root subroutine TANUPR Time average of the periodic Nusselt number, Nup U(I) S(I, J-1) V(K) *sw 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 xvi

CHAPTER I INTRODUCTION In recent years considerable attention has been given to systems containing heat exchangers operating under a wide variety of conditions. As the complexity of these systems has grown, so has the demand for more detailed knowledge of the dynamic character of their behavior. Graetz was the first to present analytical solutions for heat transfer to laminar flow in tubes. His initial work(7) considered the time-steady response of isothermal slug flow subjected to a step change in temperature at the tube wall. Graetz then investigated the same problem with fully developed Hagen-Poiseuille flow.(8) Both of these solutions are of the infinite series type and were obtained by what are now considered classical mathematical techniques. The Graetz solution for the parabolic velocity profile converges rather slowly in the entrance region so that accurate evaluation of the series solution requires a large number of the eigenvalues. These were difficult to evaluate prior to the advent of the digital computer. Leveque(l3) presented a solution for heat transfer from an isothermal flat plate to laminar flow parallel to the surface of the plate. This solution is also valid in the entrance region of a tube where the region of the fluid influenced by the wall temperature is a very thin layer adjacent to the tube wall. Sellars, Tribus and Klein,(25) using the WKB approximation,(37) have extended the Graetz solution by obtaining the larger eigenvalues. Their solution in the entrance region agrees with the Leveque solution for isothermal flat plates. Brown(2) has calculated -1

-2the first ten eigenvalues for the Graetz problem and for the companion problem of Poiseuille flow between parallel isothermal plates. Several other solutions exist for the parallel plate problem. (3,22)33) The work of Sellars, Tribus and Klein also includes the solution for a constant wall heat flux. Siegel, Sparrow and Hallman(29) present a series solution, similar in form to Graetz', for the constant wall heat flux problem which, by an integral technique, can be extended to include the case where the wall flux varies arbitrarily in the axial direction. The previously cited literature deal with the idealized problem. That is, they consider uniform boundary conditions, steady flow with fully developed velocity profiles, constant physical properties, no viscous dissipation and no secondary effects due to free convection. The literature abounds with investigations which take into consideration these effects which idealization neglects. Singh(3031) outlines an analytical method of solution for the Graetz problem including the effects of axial conduction, viscous dissipation, and energy generation within the fluid. Maslen(l5) presents analytical solutions which include the effects of compressibility, axial body forces, and variable viscosity and thermal conductivity. Several investigators have developed methods for handling the Graetz problem with a non-parabolic velocity profile. Schenk and Van Laar(23) use the Prandtl-Eyring velocity profile while Whiteman and Drake(34) and Lyche and Bird(14) use power-law profiles. Pigford(21) and Koppel and Smith(ll) consider the effect of temperature dependent density and viscosity on the velocity profile and heat transfer rate.

-3A numerical solution using the developing velocity profiles of Langhaar(12) is given by Kays.(10) The effect of non-steady or non-uniform boundary conditions has been considered by Siegel(26) who investigated the effect of arbitrary time-variations in the tube wall temperature on fully developed laminar flow in tubes. He later extended his methods to include both spatial and time dependent wall temperature (including the effect of wall heat capacity) acting on slug flow between parallel plates. (27) Siegel and Perlmutter(28) have obtained solutions for pulsating laminar flow between parallel plates with isothermal or constant heat flux boundary conditions. They have later included the simultaneous effect of a sinusoidal wall temperature or heat flux.(20) The work presented here is a theoretical and experimental study of laminar flow in a tube subjected to sinusoidally time-dependent energy generation in the tube wall. An approximate solution is obtained by the numerical solution of a set of finite difference equations derived from the energy equations for the fluid and the wall. Both the fluid and the wall equations consider radial as well as axial distribution of temperature. The experimental work consists of temperature measurement in a system which approximates the analytical model. Good agreement between the numerical and experimental results is obtained for large laminar Reynolds numbers and small energy generation rates.

CHAPTER II ANALYSIS Ao Introduction The problem being considered is the steady-periodic (time) behavior of the temperature of the system shown in Figure 1. The system consists of a constant diameter tube through which a fluid flows; the velocity profile being steady, fully developed and laminaro Heat is generated in the tube wall at a constant volumetric rate which varies sinusoidally with time. The outer periphery of the tube is adiabatic, causing all of the energy generated in the wall to either increase the wall temperature or be transferred by conduction to the fluid at the fluid-wall interface. Since the energy generation is periodic in time, the nature of the wall and fluid temperature response is also periodic; both being functions of axial distance and radius as well as time. In addition, the fluid velocity is a function of radius. In the formulation of the problem the following assumptions are madeo (a) The fluid and wall physical properties are independent of temperature. This assumption leads to some error since the fluid viscosity can be quite sensitive to temperature change. The fluid adjacent to the heated wall undergoes a significant viscosity change from the fluid viscosity near the centerline; resulting in a skewing of the velocity profile. This distortion is slight in the entrance region where the fluid temperature

q =qo( I + Sinwt ) INTERNAL HEAT GENERATION \\ x TEMPERATURE r*2 FLUID U = 2Um(I —- TEM PERATURE Ri Tf(x Sr t ) /t~~~~~ FLU i DFLUID VELOCITY FLUID WALL Figure 1. Schematic Diagram of Physical Model.

-6is practically unchanged from the uniform entrance temperature ~ except for a thin layer adjacent to the wall where the velocity is essentially zero. Downstream, where the temperature profile is completely developed, the velocity profile can become considerably distorted due to large radial viscosity gradients. (b) Axial conduction of heat in the fluid and wall are negligible. This is valid where the Peclet number (Re Pr) is much greater than 100;(24) the Peclet number representing the ratio of energy transport by enthalpy flux to that by heat conduction. For the case where the fluid is a liquid metal (low Pr) this assumption would lead to gross errors in the analytical solution. (c) The fluid is incompressible. (d) The fluid velocity is laminar and steady. This is valid when Re < 2300 and where there is no time variation of the axial pressure gradient in the fluid. (e) The problem is symmetric about the centerline of the tube. For large laminar Reynolds numbers this is a valid assumption. As the Reynolds number is reduced,body forces may become significant as compared to pressure forces, resulting in a secondary flow due to free convection. This secondary flow, shown schematically in Figure 2, causes the temperature at the top of the fluid to be higher than its predicted value while the temperature at the bottom

-7Figure 2. Free Convection in a Fluid in a Horizontal Heated Tube. Figure ect of Free Convection on a Nor-////// Figure 3- Effect of Free Convection on a Normal Laminar Temperature Profile in the Entrance Region.

-8of the tube is lower than expected. This free convection effect, caused by thermally induced density gradients, should be most evident in the entrance region where the largest radial temperature gradients exist in the fluid, The radial temperature distortion is shown schematically in Figure 3. B. Formulation The first law of thermodynamics, the conservation of mass and the conservation of momentum laws applied to the system shown in Figure 1 and consistent with assumptions a-e give the following partial differential equations for the temperature in the wall and the fluid (Appendix A gives the details of the derivations) Wall (I/r RtI I" &4 = 6 sits) Ae - ar r (II-1) Fluid (0o r~ /).E, + D r r r (11-2) which must satisfy the conditions a _ /, 9o) = 0 (symmetry) air (, a), j ) 0 (adiabatic wall)

-9/' (o, Qr) O (entrance) (~ $ /) = At (~~/9) (wall-fluid interface) The Equations (I-1i) and (II-2) and the boundary conditions being linear, the solution for ~ may be assumed to be the sum of two independent functions; one steady and the other periodic. Thus Wall P(',.r) = tS(%r) + cer( (II-3) Fluid ~t; (1} 9r) = if (3,r) - any (Xv,yr) (II-4) with the sum of *s and *p giving the exact solution for Vo Substitution of Equations (II-3) and (II-4) into (II-1) and (II-2), respectively, and separation of the steady and periodic parts yields Wall __0 - / @ _i __ = 0 ( II-5) arI- r i (6

-10Fluid (i-r~) - - =_ +~..2, ( )a7( A ar r Or (II-7) r% ara D 4 (/-r2) _=, OF (II-8) subject to the conditions a IQIS~V (N A) _ _ 9 (U,;o, s O (II-9a) a ('s L (, o) _ _ (a,> 0, ) (II-9b) (/f (o, r) =- ~- bf (o,r, ) = O (II-9d) YPsUr Di = puU ( O., =, O (II-ge) The transient-periodic solution of Equations (II-6) and (ii-8) will not be attempted here. The solution being sought is the asymptotic or steady-periodic solution, valid for e -Ho. While Equations (II-6) and (II-8) are independent of any initial conditions in the asymptotic state, nevertheless, some condition related to time must be specified to satisfy the time derivatives. Since the energy generation rate is sinusoidal, qo(l + e Sin rI), it is necessary to assume the functions

-11tpf and *pw will also be sinusoidal. This condition is imposed by assuming that these functions are of the form pe = A (,r) Sinre 4- B(%,r) Cosr (II-lo) and pw = C (,r) Slnr r - D (, r) C os r (II-ll) These equations can be rearranged into the more convenient form S = \lC/ S+ o$n 6 - where A )-: a-' -C Recombining *pf with srf and Ipw with *sw, *s = # + A B4Finally, these two expressions for temperature response can be written in the same form as the energy generation rate. That is, VVt'= R ( / t A Sl in (D&-#t)) (II-12) Vj$= s (I +. f A/?r Sn g_ (11-13)

-12where the amplitude ratios are defined by ARs =VC i D2 f The functions ARf, ARw, of and Ow are all functions of x and ro This form of the solution is most desirable when comparing the analytical results with the reduced data obtained from the experiment. In addition, the response of the wall and fluid temperature to changes in the parameters rF A and A are more easily interpreted when the solutions are in the form of Equations (II-12) and (II-13)o Co Finite Difference Equations Solution of the fluid differential Equations (II-7) and (II-8) in conjunction with a modified form (radially lumped) of the wall Equations (II-5) and (II-6) was first attempted using a variety of mathematical techniques, among which were laplace transformations, integral methods, perturbation, and the method of characteristics. Due to the complexity which arose from the time-dependent fluid boundary condition (which came from the wall equation) these mathematical techniques proved inadequate and it was decided to attempt a finite difference solution. Equations (II-5, -6, -7 and -8) are solved by iteration on the IBM 7090 digital computer after they are written in finite difference form. The steady Equations, (II-5) and (II-7), can be transformed directly into difference equations but Equations (II-6) and (It8) require

-13the application of Equations (II-11) and (II-O), respectively, before they are written in difference form. For example, Equations (II-8) and (II-10) combine to give Z)(oA /8 C OsT - L 8 S i,~ z_ + (/ - r S1 TO + "gCOST' 3AT'CQ 7O 8rs re (, r)f/s ~ @B9~ Cosos9r, arz +-r ar az r r ar which must be valid for all @0( >> 0). Consequently, the sum of all terms with Sin rG coefficients can be equated to zero and the same for the sum of all terms with Cos Ip coefficients, This leads to the two equations -Ao +./-r = - + r (11-14) D ax ae A 4 DA(II-15) Using the same procedure, Equations (II-6) and (II-11) combine to yield rC =a r r r,) (11-16) A-D = ( +a r ar? +E6 (II-17) Where before Equations (II-6) and (I-8) were coupled at the fluid-wall interface, now Equations (II-14) and (II-16) (Cos iG equations) are

-14coupled there, as are Equations (II-15) and (II-17). The derivations of the finite difference representation of Equations (II-5, -7, -14, -15y -16 and -17) appear in Appendix B; the results are shown in Table I. The system of difference equations is derived from a grid network wherein the regions (O < r < 1) (fluid) and (1 < r < R-) (wall) are divided = Ri into elements, both radially and axially, as shown in Figure 4, and the derivatives at any point are written as functions of the variables at the surrounding points in the usual manner for a finite difference problem. The equations are designed for a marching solution in the axial direction in much the same manner as in the time dimension in a transient finite difference problem. Starting at the entrance (x = O) the axial dimension is incremented while at each successive axial distance the difference equations are solved by iterating through the entire list of equations until the variables converge to within a certain per cent of their true value. Were this marching solution not used, convergence of the solution would not be possible by iteration. The Gauss-Seidel method is employed in the iteration process. With this method., during any one iteration, when the new value of any variable has been computed it is immediately used wherever that variable appears in another difference equation. This procedure is in contrast to the normal iteration process where the new values of all the variables are not used until the next iteration. The Gauss-Seidel method usually converges more rapidly than normal iteration in the regions where both methods are stable. (5) Do Computer Solutions The computer program for solving the set of equations in Table I appears in Appendix Co Also appearing in Appendix C is another

-15TABLE I DIFFERENCE EQUATIONS - LARGE r Computer Grid Network Eqn. Function Notation Location Finite Difference EQuation a. *8f(OJ) S(O) Centerline S(O) = 4N2S(1 ) + P/L.U(O) 4N2 + P/L b. A(O,J) A(O) Centerline A(O) 4N2 B(l) -(4 + P/L) B(O) + P F(o) r r Lr c. B(O,J) B(O) Centerline B(O) = - 4 N2 A(l) + + P/L) A(O) - P E(O)' r Lr General Fluid N2(1 + 1 )S(I+1) + N2(1 - 1)S(I-1) + (1 - 2)U(I) d. *sf(I,J) S(I) Interior S(I) = 21 21 N2 Point 2N2 + 1 - 2 L N2+ General Fluid 2 2 12 e. A)Ii) pt(1- 2I 2(I- - [-~e (IJ) A(I) Interior A(I)= 1 + 1)B(I+) + 1 )B(I-) - B(I) + (1 - F(I Point General Fluid N2 2 P 12 P 12 f. B(IJ) B(I) interior B(I) = - (1 A(I+) )A(I-1) + [2N + ( ]A( Poin --- 21 r 21 I L A( -1 r P Point rsf(NJ) S(N) Fluid-Wall S(N) V(O) + V(1) + (N - A)S(N-l) + AM 1 1 *sw(O,J) V(O) Interface + ) + (N - 2) A(N,J) A(N) Fluid-Wall + 1)D(l) - [A(M + 1) + N - 1)D(o) + (N B(N-l) h. A(N) C(O) - LT + C(O,J) C(O) Iriterface r( + 1) 2MP 2N B(N,J) B(N) Fluid-Wall B(N) =+ D )C() (2 + [ 2) + N - 2]C(O) - (N - 2)A(N-1) 2M n(OJ) D(O) Interface r( 8 + ) J. *sw(KJ) V(K) GeInteriora V(K) (1 +V( K+1) + (2( 2 )V(K-) + 2 Point =2 2__ 4(M+K) 2A2 General Wall A M2 2 A M2 M2 2AM2 k. C(K,J) C(K) Interior C(K) = 1(A-1 + )D(K+l) + (- M )D(K-1) 2- D(K) Wall M22M+K8) I 82 28(M+K) rs2 General Wall -AM2 _.+ M__)(, ) A+M2 M2 +2AM2CK N' 1. D(K,J) D(K) Interior D(K) = + C(K+) ( )C(K-1) + C(K) Wall r 87 28(M+K8) rI' 2 28(M+K8) - r (M) Outer Wall (1+V=( ). sw(MJ) V(M) V(M) = v(M1)+ + Perimeter (- + M - -) 8 2 Outer Wall A(M + M - )D(M) + A(M +M 1)(M-1) C(M,J) C(M) C(M) = 8 2 8aK p. D(M,J) D(M) D(M) = + 8Perimeter r( Perimeter P6( 1+5)

(Mt?.r Q~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~,J J+IO RO~~~~~ 0 O c Figur I —, Cocnt r0 etok

program for solving these equations in a rearranged form. The necessity for this second program cannot be fully understood until iteration stability is discussed later in this section, after which the remaining program is developed. Both of the programs presented are written in the Michigan Algorithm Decoder language for solution on the IBM 7090 digital computer. They are basically divided into seven parts which are (A) Input (B) Stability tests (C) Difference equation iteration (D) Convergence tests (E) Periodic amplitude and phase lag computation (F) Nusselt number computation (G) Output If brief, data are read into the computer in the form of the dimensionless parameters which govern the difference equations. Stability tests are made to insure that the solution will converge upon iteration of the equations with these data. If convergence is not assured the data are rejected and new data are read in. If convergence is assured, the program proceeds to compute a set of parameter-dependent constants which frequently appear in the difference equations. Iteration of the equations then proceeds, starting from an initial estimate that all the variables are zero, at the first axial location Ax downstream from the entrance. At the end of each iteration, tests are made to check whether all the variables at this axial location have converged to within a preset per cent of their true value. If this is not true the program returns

to the iteration of the difference equations. Provision is made to limit the amount of iteration at a particular axial location to eliminate data which provide a stable but slowly convergent solution. If all the convergence requirements are met the program jumps out of the iteration loop and proceeds to compute the periodic amplitude and phase lag of the periodic solution at each point of the grid network for the current value of x, The local Nusselt number is computed next, after which all of the results are printed out. Then the program increments x to the next axial location downstream where the entire iteration process is repeated. When the solution has marched down the tube to a specified length L the solution is terminated and new data are read ino The process continues until the input data are exhausted. A few features of the program require some detailed explanation. 1. Stability Tests The difference equations are a set of simultaneous linear algebraic equations whose coefficients form a 3(M+N+1) by 3(M+N+l1) determinant; M being the number of grid elements radially in the wall (not including the wall-fluid interface) and N the number of grid elements radially in the fluid (not including the centerline element). Each row of the determinant represents the coefficients of one of the difference equations while the diagonal elements are the coefficients of the dependent variableso When the grid network is finely divided, i.e., when M and N are large, a majority of the elements in any row of the determinant are zero, making the Gauss-Seidel iteration method well suited for solving the set of equations.

-19A sufficient condition for iteration stability is that the absolute value of the element on the main-diagonal of the coefficient determinant be greater than or equal to the sum of the absolute values of all the other elements in the same row.(5) That is, for the i-th row. 3(M+Nl/) |a,, 66| 2 E | Q (II-18) As an example of the iteration stability to be derived by satisfying Equation (II-18), consider the simple set of equations which are satisfied by x = 1, y = 2. The determinant of the coefficients is 2I which satisfies Equation (II-18). Solving for the variable whose coefficient is the main-diagonal element?c = 4-;2,

-20Any initial values can be employed. Assuming these to be 0 and iterating: c = Z-o =z S-2= 3h ~Z 2 X_ 4 -z 4 /- /3 - ZC S Upon further iteration closer approximation of the true solution would be achieved. If, on the other hand, the rows of the determinant were switched so that the determinant was Equation (II-18) would not be satisfied. Again solving for the variable whose coefficient is the main-diagonal element AZ= 4?4-% and iterating from zero initial values ~=4-0=4 %- 5-8 =-3 _- -I = /O _ 4 3= - 34 -- tOo

-21The solution is now diverging. Equation (II-18) is equivalent to saying that the absolute value of the coefficient of the dependent variable in any difference equation must be greater than or equal to the sum of the absolute values of the coefficients of the independent variables in that equation. Since these stability inequalities can be determined directly from the equations in Table I, the entire determinant of the coefficients need not be written. The stability inequalities are listed in Table II. All of the inequalities in Table II need not be tested before iteration of the equations. Inequalities a, c, e, g and i are always satisfied so that convergence is assured for the steady equations regardless of the grid network dimensions. Inequalities b, d, f, h and j depend upon F, 6, B and A as well as M, N and P so that stability is confined to certain regions as defined by the inequalities. Some of the inequalities from the periodic equations are more strict than others and, consequently, by satisfying inequalities b and f, all the remaining inequalities are satisfied. It is these inequalities which are tested in part B of the computer program. Starting with large values of N and M, these parameters are reduced until the inequalities are satisfied. Lower limits are placed on N and M such that if stability has not yet been achieved the input data are rejected and a statement is printed telling whether instability is caused by the wall or fluid equations. The lower limits on N and M are necessary because a very course network in either the wall or fluid gives rise to a sizeable error in the finite difference expression of the derivatives. This is discussed in detail in Appendix Bo

-22TABLE II STABILITY INEQUALITIES - LARGE r EQUATIONS Equation and Inequality Location Stability Requirement ao Steady- 4N2 Centerline 1 >4N2 + P/L b. Periodic- > 4N2 + 4N2 + P/L Centerline r r c. Steady- I General 1> 2N2 L N2 Point-Fluid 2N2 + (1 - 2 + l do Periodic- 12 General >2N2 2N2 + (_) Interior - 2N + p Point-Fluid eo Steady-Wall- 1> + N1 + Fluid Interface = AM + i) + N - + ) + N2AM + }} + 2( f. Periodic-Wall- 1 (N 1) Fluid Interface r 2( M 2N go Steady-General 2M2 62 Interior 1 >2M2/2 Point -Wall ho Periodic-General 4 Interior 1 - Point -Wall io Steady-Outer M/5 + M - 1/2 Perimeter Wall M/16 + M - 1/2 j. oPeriodic-Outer 2A(M/6 + M - 1 Perimeter Wall 2MM

-23In the fluid equations, stability may not be achieved, even when N is reduced to its lower limit, when r is small. This would seem to limit the solution to large values of Fo However, inspection of Equations (b, c, e and f).in Table I shows that rearrangement is possible, resulting in a marked shift in the stability requirements for these equations. For example, taking Equation (e) and solving for B(I) /B(= Nz2(/Si)lr 8(r. NS(/-1)BI(-) TA ( %('-N)F) (II-19) The general stability inequality, Equation (11-18), now requires that Fo!h otact ca w For the most accute case, when I = N-i, L 2/J - N ( I 2r (II-20) Equations (b, c and f) in Table I must be rearranged similar to Equation (II-19)o Now with a given minimum N (minimum from the standpoint of difference derivative error) and r small, P can be increased until Equation (II 20) is satisfied. All of the rearranged fluid equations will now be stable. Mathematically, what this equation rearrangement amounts to is an exchange of rows in the periodic fluid section of the determinant of the difference equations coefficients to provide a main-diagonal

element of sufficient size to satisfy Equation (II-18)o The one disadvantage of this procedure is that every time the number of axial elements P, is increased to satisfy Equation (II-20), more increments in Ax, P, are necessary to march down the tube to a given x, JL P In the wall, when both r and the dimensionless wall thickness 6 are small, the wall periodic stability inequalities cannot be satisfied. It is not possible to shift rows in the coefficient determinant to give a suitable main-diagonal element as shown by an inspection of Equations (k, 1, n and p) in Table Io What is done is to assume that the entire wall can be lumped radially, oeo,, that there are little or no variations in wall temperature in the radial direction at any x. This approach requires a reformulation of the difference equations where the last half element in the fluid and the entire radially-lumped wall are considered to be at the same temperatureo This appears at the end of Appendix Bo The second computer program in Appendix C is a program which combines the rearranged fluid equations similar to Equation (II-19) and the lumped wall approach. The difference equations are listed in Table IIIo Satisfying Equation (II-20) insures stability for the fluid equations. Applying the general stability inequality to the wall equation gives - ( &ei+\7) ~ N-Z (II-21)

-25TABLE III DIFFERENCE EQUATIONS - SMALL r Computer Grid Network Eqn. Function Notation Location Finite Difference Equation a. *sf(OJ) S(O) Centerline S(0) = 4s (1) + P/L U() 412 + P/L b. A(O,J) A(O) Centerline A(O) = 4N2A(l) + rB(0) + P/L E(O) 4N2 + P/L c. B(0,J) B(O) Centerline B(O) = 4N2B(i) - rA(O) + P/L F(O) 4N1 + P/L General Fluid N2(l + I)S(I+) + N2(l - 1)S(I-1) + P I2 sd. *f(I,J) S(I) Interior S(I) = 2121 L I Point 2N2 + (1 - ) General Fluid N2(1 + 1)A(I+l) + N2(1 - 1 )A(I-1) + rB(I) + (1 - 2)E(I) e. A( I,J) A(I) Interior A(I) = 21 21 L Point 12 Point 2N2 + l - General Fluid N2(l+ = I)B(I+l) + N2(1 - 2)B(I- 1) - rA(I) + (l - )FI) f. B(I,J) B(I) Interior B(I) = 2 2 Point 2N2 + L(1_ )2 Fluid-Wall g. *sf(NJ) S(N) Interface- S(N) = S(N-I) + 82 + 28 Lumped Wall 2N-1 Fluid-Wall h. A(N,J) A(N) Interface- A(N) = (2N-1)B(N) + (2N-1)B(N-1) Lumped Wall r( + 28 + 1 B(NJ) B(N) terface- B(N) (2N-1)A(N) - (2N-1)A(N-1) - (52 + 28)e i. B(N,J) B(N) Interface- B(N)2 + 1) Lumped Wall r(82 + 2 N+

-26In part B of the second program, Equations (II-20) and (II-21) are tested. Starting with an initially large value of N, N is reduced until Equation (II-21) is satisfied. Then Equation (II-20) is tested and P increased until the inequality is satisfied. An upper limit is placed on P so that tx,, will not become too small and require excessive computer time to march to x = L. When 5 becomes very small, the difference between the minimum r to satisfy the first program and the maximum value. of r to satisfy the second program becomes appreciable. For this case a third program was used. Here the lumped wall portion of the second program was employed while the fluid equations were those from the first program. This program is not presented here since it is identical to.the second program except for the fluid equations and fluid equation stability inequality from the first program. For larger wall thicknesses the regions of convergence for the first and second programs overlapped and the third program was not necessary. 2. Convergence Tests After iterating through the set of fluid and wall difference equations a test must be made to check whether the solution has converged to an arbitrary degree. This is accomplished by testing whether each of the variables [A(I), S(I), V(K), etc.] has changed a preset per cent between one iteration and the next. After an iteration has been completed, the current values of all the variables are stored for this convergence test. For example, if convergence has not been achieved, A(I) is stored in G(I). Iteration then proceeds again through the entire

-27set of difference equationso Now A(I) has a new value and a test is made to see whether A (r)- (r) - ALPHA (II-22) A(I) The smaller ALPHA is, the closer one comes to the true value of A(I) but more iterations will be required to satisfy Equation (II-22), When this convergence inequality is satisfied for all the variables the iteration process ceases. In general, for a given ALPHA, the closer the stability inequality comes to not being met, the more iterations will be required to satisfy Equation (II-22). Provision is made in the program to count the number of iterations that have taken place and when this exceeds a preset number, iteration is terminated, a statement is printed telling where the slow convergence exists, and new data are read in. This also acts as a backstop for the stability tests in the improbable event that some instability-producing data might pass through the stability tests to the iteration stage of the program. In any given iteration method for solving simultaneous equations the closer the first guess is to the true values of the variables the faster will be the convergence to the true valueso In the program the initial iterations take place at the first axial location downstream from the entrance, where all the variables are zero. Therefore, as a first guess for this location all the variables are equated to zero. For each succeeding axial location the initial estimate is the final

-28value of each variable at the previous axial location. This is accomplished by letting each of,the variables [A(I), S(I), V(K) etc.] retain their final values at the end of the iteration process at the previous axial location. When x is incremented and iteration is restarted the final values of the variables at x - Ax are in memory and thus become the first estimate at x. The axial derivatives of the variables, e.g. in the wall and in the fluid near the wall are the greatest at the entrance and they decrease as x increases. Thus, at the first axial location the first guess of zero for all the variables is much too low and convergence may require hundreds of iterations. As x increases this drops off and after three or four increments in x the required number of iterations may fall to as low as five. 3. Nusselt Number Computation When the convergence inequality, Equation (II-22), has been satisfied for all the variables at a given axial location the local Nusselt number is computed for that x. The local Nusselt number is defined by heA(7Ere - Tr Tmm) A kfA (II-25) where T.,, the mixed-mean or bulk temperature of the fluid, is defined by 4'>4/OlcUrnab (/I )7- rr* r"

-29which becomes, for the case where all physical properties are considered constant r* Trr = ~' (`-eV 1)r 9dr In terms of the dimensionless temperature and the dimensionless variables this is i j ( -r dr or 7-nn = _ _ + 1 o (I 25) Substituting Equation (II-25) into (II-23) yields..p P' + kf'(-r )d / 4 \;r rI Converting all temperatures to their dimensionless form h (._Z' r (I-r )dr' a hp j~~~~rr/ r~ ~ -

-30or or' r=i NuP = -r r (11-26) r' (/- r )dr Equation (11-26) can now be expressed in terms of the computer solution where \d t'= I - - = S(N) - S(N-I) 4_ A (N)- A( 1 (II-27) Nrel, SI) - A(N) Sn 4-9 ) 8(N) Cos(1 (II-28) it oaf' -' )dr = R 1? (So) 4- A (O) S'n I'& + B(O) CO( S S t 4' (svl- ) r/.17

-31The first term of Equation (II-29) represents the contribution of the centerline element to the mixed mean temperature and the summation represents the contribution of the remainder of the elements. The summation ceases at I = N - 1 since U = 0 at r = 1 Employing Equations (II-27, -28, and -29) in (II-26) yields IV[s r(N)-S(N-) A (sN)r)A( /V-/) Si T 8B(0 - 8I CosAN Nuo = L_ / _ I 5 ( N)-.(8% 4- 4 ) ( N —) ZN1 N E 4I),I -2 A( |. +- 8 I V -/3 )O) st_ -, - (II-30) The steady Nusselt number (the Nusselt number which results from a steady volumetric generation rate qo) is S (N) - S S N-l)) This also is the steady part of Nu Inspection of the set of fluid periodic equations, Equations (II-14) and (II-15), shows Lim A (I) = L', B(r) =O (-= / O.. A) Thus, from Equation (II-30) Li-n Nu - NuA E-* o

-32which was expected from Equation (II-31). The time average of Nup is Nus; that is 17r2 t+ r,/ N,,p de~Nu dN This is readily verified by integrating Equation (II-23) between the limits @ = Ty T + (w1 ( T*R, T Tm) de =kf / ( which, in terms of the dimensionless variables, is j 2,gf (r rr= _-r Nu, & (/- _rl) dr Since 4f = Isf + A Sin rG + B Cos rG the time average of *f will be \rsf Thus _Nup =,_ r~l - S Ql ( r r-32)

-33or \Up hN ("-33) This is to be expected since the time-average of the response of a linear system to a periodic input about a mean value is the response of the system to the mean value. Numerical integration of Nup from 9 = 0 to Q = 2?T in t/6 increments was included in part F of the program for several check runs on the computer. It was found that Equation (II-33) is satisfied within the accuracy of the numerical integration.

CHAPTER III NUMERICAL RESULTS The programs were run on the IBM 7090 digital computer at the University of Michigan Computing Centero Runs were made varying one parameter at a time and marching to a maximum axial distance of o012. This relatively small limit on x was decided upon in view of the fact that little further knowledge would be gained from the results at large values of x and, more important, there was the matter of economicso With computer time costing in the neighborhood of $8O00/minute, the most practical approach was to limit the size of x and run the programs with a wide variety of system parameters while still keeping the computer bill at a modest level. A majority of the solutions were run with the temperature dependent fluid and wall parameters, P and A, set for the system of stainless steel-water (the materials used in the experiment) at 60"F. The dimensionless frequency, r, was varied from 50 to 6000 with dimensionless wall thickness of from 0o2 to lo0o A few runs were made varying $ and A to investigate the effect of materials on the solutiono The results presented here are not intended to be all-inclusiveo Such a coverage would be ponderous and would require excessive computer time. Rather, the intent here is to present the general characteristics of the solution and to point out trends which appear to be functions of the problem's parameters. The results are presented in two partso steady and steady-periodic. -341

-35Ao Steady Solution The steady solution is the response of the system to a uniform steady volumetric energy generation rate qo in the tube wall. The steady dimensionless temperature, 4s, is a function of x, r, 8 and P/Ae Figure 5 gives the steady Nusselt number, defined in Equation (II-31), as a function of axial distance. These results compare favorably with Sellars, Tribus and Klein(25) whose solution is valid for a constant wall heat flux, which is what the steady solution presented here reduces to. Sellars, Tribus and Klein give three expressions for the Nusselt number: a series solution valid for all x but slowly convergent in the entrance region, a solution valid for small x and the asymptotic value which is 4.36. The solution valid for small x, which is Nus = 1.6393 X is used in Figure 5. Also shown in Figure 5 is a correlation given by Kays(lO) for his numerical solution employing a fully developed laminar velocity profile. This correlation is expressed as Nus = 4.36 - 0.04-W/% / 4- 0.00247/ The time-average of the periodic Nusselt number, Equation (II-30), was taken at large and small frequencies and compares excellently with the steady Nusselt number. Figure 6 is a plot of the periodic Nusselt number, as a function of time, at a high frequency. As r is decreased,

4 KAYS 10 Nus==4.36 + 94602 x' 2 14.02 (251 V SELLARS,TRIBUS & KLEIN - Nu 16393 NUMERICAL 10 Nus' 8 7 6 ON 5 4 3 0 2 4 6 8a 1 -&(Re Pr) xO1 R Figure 5. Steady Nusselt Number versus Axial Distance,

I Il8=0.2 15.870 r=6000 x=O0.001 E =0.25 5 = 1. 1 15.860 A = 26.89 15.850 15. 840 q 15.830 NuNup 15.820 Nu, =15.84535 Nup= 15.84532 0 2 t3 wt. Figure 6. Periodic Nusselt Number versus Time.

-38the amplitude of Nup approaches e; that is r70N Ao X + (UN),ntJ while the phase angle of Nup approaches 0 in the limit. Figure 7 gives the variation of the steady temperature at the fluid-wall interface with axial distance for several wall thicknesses. The curves are somewhat misleading in that an increase in wall thickness appears to produce an increase in temperature at the fluid-wall interface. This is true if one considers the problem from the standpoint of a constant volumetric generation rate qo; the greater the wall thickness, the greater will be the total energy generation rate in the tube and the higher will be all temperatures in the tube and fluid. However, considering the case where the total generation rate is constant for all wall thicknesses, the energy generation rate/unit volume of tube decreases as wall thickness increases and the product 25 (Sro23 ~ )) &T- To 2 / is a constant. This simply means that the actual temperature difference between the entrance and any dimensionless location in the system is independent of wall thickness when the total energy generation rate in the tube is constant. Figure 8 shows the distribution of fluid temperature with radius for various values of axial distance. The temperature profile

-393.6 r =I.0 3.2 043 2.8 2.4 2.0 1.6 0 2 4 6 8 I 0 12 Ri Figure 7. Steady Temperature at r = 1.0 versus Axial Distance for Various Wall Thicknesses.

1.0 0.9 = 0.4 A=0. 043 0.8 r 0* 7 Ri X 0. 001 0.6 ~~~~~~=0.004 0.6~~~~~~~~~~~~~~~~~~~~~~:0.004 0 0. 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 +I'( I) Figure 8. Steady Fluid Temperature versus Radius for Various Axial Distances; The Development of the Steady Temperature Profile in the Fluid.

-41can be seen developing as x increases. This development continues up to x = 0.25 after which all temperatures in the fluid and wall increase linearly with x and the temperature profile remains unchanged.(25) This is the point at which the steady Nusselt number reaches its asymptotic value of 4.36. Figure 9 shows the development of the temperature profile in the wall, comparable to the developing profile in the fluid in Figure 8. The actual wall temperature difference, remains relatively constant as x increases. In Figure 9 the "normalized" wall temperature difference appears to decrease as x increases. This is due to the fact that the normalizing factor, tsw(r=l), increases as x increases. Figure 10 gives the steady temperature at the fluid-wall interface as a function of axial distance for various values of the ratio n/A. A plot of the normalized steady temperature as a function of radius for various P/A, similar to Figure 8, shows that the development of the temperature profile in the fluid is independent of P/A. Figure 11, on the other hand, shows the temperature profile in the wall changing markedly with P/A. In addition, with 4sw(r=l) decreasing as P/A decreases, the actual wall temperature difference decreases rapidly as P/A decreases. Thus, when the tube material has a large thermal conductivity ( such as copper ) and the fluid thermal conductivity is moderate to small, the temperature drop across the wall will be slight and the assumption that the wall may be lumped radially will produce an accurate solution.

1.4 8=0.4 m:0.043 g3 A 1.3 1.2 R I. I.0 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08''s(.) Figure 9. Steady Wall Temperature versus Radius for Various Axial Distances; The Development of the Steady Temperature Profile in the Wall.

-431.2 i, 0.043 1.0 0 0.6'.0 0.4 r —-.0 o 06 0.2 0 2 4 6 8 10 12 x (Re Pr )"xlOI Ri Figure 10. Steady Temperature at r = 1.0 versus Axial Distance for Various 5/A.

I I 1.4 1.3.0 A o. o43 A r 1.20 A -0. 010 1.1 x a0.012A 1.0 I I.i' —~ ~.P~ 8:0.4 ~ ~ ~ 0.02 1.0LO~~~~~~~~~, 1.01 1.020 Fge Seyleer.arvsuRdis B~o.4 1.0~~~~~~~~~~~~~~~ I I 1. i 1.01 1.02 Figure 11. Steady Wall Temperature. versus Radius for Various B/A.

-45B. Steady-Periodic Solution The steady-periodic solution, a function of r, x, E, r, 6, f and A, is the response of the system to an oscillating energy generation rate in the wall, qoE Sin r@. Inspection of Equations (II-14) and (II-16), which are coupled at the fluid-wall interface, shows that as r approaches zero, B and D approach zero since there is no energy generation term in either equation and the derivatives in the radial direction are zero at the centerline and at the outer perimeter of the tube, Equations (II-15) and (II-17), on the other hand, are influenced by the energy generation term Pe in the wall equation and consequently their solutions are not trivial in the limit. Comparison of Equation (II-15) with Equation (II-7) and (II-17) with (II-5) shows that thus Lamw A A = s.im\/ArB - A = 6'Sf (111-2) Lim\lc C D' - C - 6 5 ran~~~~~~~~~~~~~~~zze

-46Further Lm -L ro r A L-m ot = L nrn (~~"n A) ~- c (III-3).r — W- O rho — ~ O The steady-periodic results are presented in the form of phase angle, 0, and amplitude ratio, AR = 1A2 + B2/Esf. The periodic amplitude, IA2 + B2, is directly proportional to E as Equations (II-14, -15, -16, and -17) show. As a consequence, the amplitude ratio is independent of E o Figure 12 presents the temperature phase angle at the fluidwall interface as a function of axial distance for various values of wall thickness. The curves show that the wall exerts a significant influence on the temperature there; generally, the smaller the wall thickness the smaller will be the temperature phase angle in the wall and at the fluidwall interface. The phase angle shows a tendency to oscillate with x in a slightly underdamped manner with an amplitude which decreases as 6 increases and at a frequency independent of 6. The degree of damping appears to be dependent upon 6; increasing as 8 increases. Figures 13 and 14 give the phase angle in the fluid at r = 0.9 and r = 0.8 respectively. These curves show the same tendencies as those in Figure 12 with the maximum phase angle occurring at a larger x as r decreases and with less damping at the smaller radii. At the same time, the period of oscillation in the x direction increased as r decreases.

-471.55: -~:o.6~8=06 1.501.0.4 r= l.lO -=1.16 1.45- &= 26.9 1.4 I I I 0 2 4 6.8 10 12 * -1 3 - (RePr) x 10 Figure 12. Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r = 1.O.

-482.5 2.4 2.3 9 1.0 2. -8 0.8 c_ r-0.9 a. 0.4 r.200o_ 0-o.2 B~1.1 Ax.26.9 2.0 1.9 0 2 4 6 8 10 12 l-(RePr) xIO Figure 13. Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r a 0.9.

-493,6 3.4 3.2 E =1.0 3.0 -.| z r=200 4 /// \ 8=0.4 r =0.8 4l.16 2.8 8=0.2 re~~~~~~~~~~~~A 2 =26.9 2.6 2.4 2.2 2.0 I I I I I 0 2 4 6 8 10 12 R (RePr)x 103 Figure 14. Fluid Temperature Phase Angle versus Axial Distance for Various Wall Thicknesses, r = 0.8.

-50Figure 15 gives the temperature amplitude ratio at the fluid. wall interface as a function of axial distance for various wall thicknesses, comparable to the phase angles in Figure 12. The function IA2 + B2 is zero at the entrance but increases rapidly to a steady value invarient with x (for all but very small r) just inside the entrance of the tube. As a consequence, the amplitude ratio varies with x as l/*s and, with *s increasing linearly with x for x > 0.252 the amplitude ratio approaches zero in the limit as x becomes large. Also, at a given x the function IA2 + B2 increases only slightly with 6, unlike 4s where 2s is practically constant. Thus, AR des52 + 26 creases continuously as the wall thickness increases. At a constant wall thickness the amplitude ratio increases as x decreases and approaches a limit of 1.0 as x approaches zero. This is verified by noting LAim AR = Lim \ oA ( -4) From l'Hospital's rule Lim AR =I (A a (% A= \]=o _/'/ A. B L I(, A 8'' as J~ (III-)

'0'T = I'sassa'Ul'Tqi, T'-e snOTtA ioj aoou.sTCE ~TVTXV snsgaaA OT- T apnl4TTdu a~ n.vadmuaj, pTnT,'T5 a'n.Tf 01 x (Jd 98 ). 01 8 9 1b 0 Z'O 9'0 0'1-, I80;I // 0= 91z~~org 2~~~~0'1 \ 9'01 \ 009 -J V 01 -= O I - -.

-52Now, from Equations (11-14, -15, -16 and -17), the only energy generation term in (II-14) and (II-16) is FA and rC. Since [mA - LimOC 0 tB 6o and since both are small for small x, B, -x, D and ax are — O when x is small. Equations (II-15) and (II-17) have the generation terms rB and ID + PE. Thus A I SMALL ~a > C aJ Therefore -aA Lim AR -iyr -o With B = O, A = Ec:sx and thus Lim A?: / (III-6) for all.v At the same time Lim = LL 0 (111-7) 1 —~- 0 —,O Awhich is also valid for all 5 Figure 16 shows the variation in the temperature phase angle at the fluid-wall interface with axial distance for various frequencies. As r increases at any x the phase angle increases and approaches c/2 as a limit as r approaches o. The phase angle oscillates with x at a frequency which increases with r and is damped at a rate proportional to r. The frequency of oscillation in the x direction also

-53L54 r = 600 1.50 r=500 r = 400 r = 300 1.46 z r=200 1.42 1.38 1.34- r =1.0 E-1.0 8=0.2 /G=1.16 1.30 A:=26.9 I I I l I 0 2 4 6 8 10 12 X- (RePr)- x 103 Figure 16. Fluid Temperature Phase Angle versus Axial Distance for Various Generation Frequencies.

-54appears to decrease as x increases. This is shown in Figure 17 which is a magnification of the curve for r = 600 in Figure 16. Here the slightly under-damped oscillation is evident, with the period of oscillation increasing with x as stated above. Figure 18 gives the temperature amplitude ratio at the fluid-wall interface, corresponding to the phase angles given in Figure 16. The amplitude ratio at any x increases as r decreases and tends towards a limit of 1.0 as r goes to zero. Here, as in Figure 15, IA2 + B2 remains constant for all values of x except at the entrance. Again, with *sf continuously increasing with x, the amplitude ratio decreases with x and tends toward a limit of zero as x becomes large, except at r = 0. As x approaches zero, the amplitude ratio again approaches 1.0 as shown in Equations (III-4) to (III-7). When r is large, Figure 18 appears to show AR approaching a limit less than 1.0. This is due to the fact that the energy generation term rA in Equation (II-14) does not become small, with r large, until x is very small. Thus the tendencies outlined in Equations (III-4) to (III-7) are valid only at smaller and smaller x as r increaseso Figure 19 presents the fluid temperature phase angle at various radii as a function of axial distance, for a small r-. At and near the fluid-wall interface the effect of the wall predominates in that the phase angle appears to oscillate only slightly in the axial direction. As r decreases at a given x the phase angle increases. At a given radius, phase angle increases with x, reaches a maximum value and oscillates about a mean value in a slightly underdamped manner. As r decreases the maximum phase angle occurs at a larger axial distance and

-55r=l.0 r=1.0 8=0.2 1.510 -'a,=1,.16 A = 26.9 1.5090 1.5080 0 2 4 6 8 10 12 Rx (RePr)- x 103 Ri Figure 17. Fluid Temperature Phase Angle versus Axial Distance.

.30.25 r = 1.0:= 1.0 8= 0.2,8= 1.16 A= 26.9.20,.15 200.05 300 400 600 0 2 4 6 8 10 12 X* (RePr)lx 103 Ri Figure 18. Fluid Temperature Amplitude Ratio versus Axial Distance for Various Generation Frequencies.

'001 = J'T ITe SlnoT lo;A J ousas'a -T.xV S1 A alT.WU asGq amaauL PTn'T'*6T an'l:T O1 xl(_ dSa) x, 01 8 9 6'9Z=V 91I'l = 0'1.:. 0'O= = 68' ~ ~Z00001= 1=, JL= 9'0 86'0 / / I - l -

-58the period of oscillation decreases. In Figure 20 the fluid temperature phase angle is plotted against axial distance for various radii with a medium-range value of r o Here the maximum phase angle at any r is reached at a smaller x, as compared to Figure 19, and the steady value about which the phase angle oscillates is larger, In Figure 21 the same curves are shown for a large r. The same tendencies as noted for small and medium r occur, with very little oscillation at the larger radii. In Figure 22 the amplitude ratios corresponding to the phase angles in Figure 19 are shown. With VA2 + B2 again invariant with x, after the entrance, for all r Lrn AR =O The form of the function representing the amplitude ratio appears to fall into three classes, depending on r. For small r the numerical solution appears to vary as AR= Kce (r <0.) (iii-8) Near the fluid-wall interface the amplitude ratio has the form AR = t,( r > o.8,J % > o.0z) (III-9) In the range (0.5 < r < 0.8) a complex function makes the transition from Equation (III-8) to Equation (III-9). At a large radius (near the wall) the curves cross and recross one-another for reasons to be discussed in conjunction with Figure 24, amplitude ratio as a function of radius at various axial distances. In Figure 23 the temperature phase angle in the fluid is plotted against radius for various axial distances. At any x the phase

-593I r r-0.5 r = 500 E= 1.0 5_r 8= 0.2 I // o~~~~~~~~~~r.= 0 2. = 0.8 r= 0.9 0 I I ~ I I 0 2 4 6 8 10 12 * -I 3 Figure 20. Fluid Temperature Phase Angle versus Axial Distance for Various Radii, r - 500.

-60r= 0.5 5 I 57r 2 2 _ r= 0.7 c3 7r r- r0.9 r= 1.0 T 2 ~r 1500 E= 1.0 8 0.2 = 1.16 M= 26.9 0 I l I I 0 2 4 6 8 10 12 (Re Prl x lO Figure 21. Fluid Temperature Phase Angle versus Axial Distance for Various Radii, r X 1500.

-612.8 1.. r= 0.5 r = 0.2 2.4 r: 0.3 2.0 r=0.9 1~r=.6~~~=0 i 0.22 0.4~~~~~~~~~~~8 P0.2 0.4 1.16 = 26.9 0 I 0 2 4 6 8 10 12 R (RePr) X 103 Figure 22. Fluid Temperature Amplitude Ratio verses Axial Distance for Various Radii, r w 100.

1.0 0.8 K =.002 K =.004 x =.006 0.6 x =.008 Fr 200 Ri I E 1.0 0.4 8 0.2 A. 26.9 0.2 0 0 2 Of (RADIANS) Figure 25. Fluid Temperature Phase Angle versus Radius for Various Axial Distances.

angle increases as r decreases and reaches a maximum at r = O. Near the centerline the rate of change of Of with x or r is practically invariant with x while at the fluid-wall interface the phase angle varies little with x (as was noted in Figure 16) and the rate of change of of with r increases only slightly with x. As r increases the curves retain the same general form and characteristics, with the magnitude at any location increasing as r increases. Figure 24 presents the temperature amplitude ratio plotted against radius, corresponding to the phase angles shown in Figure 23. The general form of these curves.is governed by the relationship which exists between IA2 + B2 and tsf' In general VA"74- ASc aAR Dr ar (III-10) In Figure 24 AR >(ar foau, <A (III-11) AR at al As r defor all x. Thus, < 0 at the fluid-wall interface. As r decreases a point is reached at which AR is equal to the ratio of the derivatives and an inflection point occurs. The radius at which this takes place is a function of x. As radius decreases further, the inequality in Equation (III-11) is reversed and the amplitude ratio

1.0 rP 200 621.0 0.8 - 80.2 ~81.16 A 26.9 0.6 r Ri ~0.4 X =.010 X =.008 X =.006 0.2 X=.004 X =.002 0 - 0 0.2 0.4 0.6 0.8 1.0 1.2 ARxlO Figure 24. Fluid Temperature Amplitude Ratio versus Radius for Various Axial Distances.

increases. At the centerline another inflection point occurs due to the derivatives on the right-hand side of Equation (III-10) being zero there, due to the conditions of symmetry imposed at the centerline. It is possible to have a condition where the amplitude ratio either increases or decreases continuously from r = 1.0 to the centerline. This can be seen in Figure 22 for x =.001 and x =.003. At x =.0015 the amplitude ratio is invariant with r, indicating that IA2 + B2 is proportional to *sS for all r. The fluid temperature phase angle is shown as a function of frequency for various axial distances in Figure 25. The solution has the characteristic that the curve for x and x/2 cross at a certain frequency and intersect again at a higher frequency - the same frequency where the curve for x/4 crosses the curve for x/2. The same pattern repeats at higher and higher frequencies where the curves for x/8, x/16, etc. cross. For a larger radius the same characteristic occurs but at successively lower frequencies for each intersection and with smaller amplitude variation from one curve to the other. The temperature amplitude ratio is plotted against frequency in Figure 26. The curves exhibit a resonant characteristic which diminishes as x increases while the region in which resonance occurs is independent of axial distance. The frequencies at which the maxima and minima are found are also independent of axial distance. For larger radii the same characteristics occur with smaller oscillation amplitudes.

57 I I 2 r 0.6 Es.O0 80.2.8x1.16 371 2r _.= = J9R 1 a X =.006 2 0 100 200 300 400 500 6 wRi2 raf Figure 25. Fluid Temperature Phase Angle versus Generation Frequency for Various Axial Distances.

9 r=O.8 1.0 8=0.2 / 31.16 A =26.9.7 X =.0015.6 X 00 3 X =.006 X=.012.5 AR.3.2.$I 0 100 200 300 400 500 600 r wRi2 af Figure 26. Fluid Temperature Amplitude Ratio versus Generation Frequency for Various Axial Distances.

CHAPTER IV EXPERIMENTS A set of experiments were performed to investigate the validity of the assumptions made in the analysis and to provide a rough check of the.entire numerical solution. In particular, the experiments were designed to determine how the temperature dependent physical properties of the fluid and the free convection phenomenon effect the temperature distribution, as predicted by the computer solution, at the outer perimeter of the tube. Temperatures were measured along the outer wall of the tube by means of thermocouples spot-welded to the wall while simultaneous measurements of entrance temperature, energy generation rate in the tube wall, and fluid mass flow rate were made. The experiments were performed in the Heat Transfer and Thermodynamics Laboratory of the Department of Mechanical Engineering in the University of Michigan. Ao Apparatus The apparatus used in the experiment is shown in Figures 27 and 28 and schematically in Figure 29. It consists of an open fluid circuit which contains a constant head reservoir, a line heater, an electrically heated test section, a weigh tank, a DC power supply with remote controls, and associated instrumentation. 1o Test Section The test section was a horizontally oriented, 5/8 inch BWG 20, 304 stainless steel tube, 10 feet in length and insulated on the outside -68

~~~~~~~-6 - Figure 27. Experimental."ppara~tus Showing Test Section.

-70Figure 28. Experimental Apparatus Showing Reservoir Tower and Controls.

POWER SUPPLY VENT FLOAT G VALVE SPEED GEAR MOTOR REDUCER BOX AUTO TRANSFORMER RESERVOIR I N WATERR LINE HEATER CONTROL IMMERSION WATER HEATER SOFTENER TEST SECTION WEIGH TANK Figure 29. Schematic Diagram of Experimental Apparatus.

-72by one inch thick fiberglass insulation. The initial six feet of the test section provided a sufficient entrance length to insure fully developed laminar flow into the last four feet of the tube. This last four feet of the test section was heated electrically by a DC power supply described below. The power was delivered to the test section through two copper electrodes, spaced 40 inches apart, silver soldered to the test section, Two pressure taps were installed in the test section, one four inches upstream and the other four inches downstream from the heated section, to provide for pressure drop measurements. Both ends of the test section were electrically insulated from the remainder of the fluid circuit by short sections of rubber tubing. 2. Fluid Supply System The fluid used in the experiment was water as received from the Ann Arbor City mains. The water was passed through a softener prior to entering the reservoir which was located on a tower above the test section. The reservoir contained an inlet control valve which maintained the free surface of the fluid ten feet above the test section. Atmospheric pressure was maintained at the free surface of the water by venting the top of the reservoir. The fluid flowed from the bottom of the reservoir to the remainder of the circuit by gravity. 30 Line Heater A General Electric 10OKW, 230VAC immersion line heater was installed in the circuit at the exit of the reservoir. Control of the power delivered to the heater was provided by two General Radio 115V

-73variacs with a continuously variable output from 0 to 1000 watts. The heater was used to control the fluid temperature at the entrance of the test section and to provide a variety of test section temperatures for calibration of the thermocouple instrumentation. A mixing baffle, shown in Figure 30, was located at the exit of the heater to insure a uniform heater exit temperature. This temperature was measured by a thermocouple probe inserted along the centerline of the circuit tubing immediately downstream from the mixing baffle. 4. Power Supply Electrical power was dissipated in the tube wall to provide the desired heat generation. The power was supplied by a transformerrectifier unit manufactured by Hansen-Van Winkle Co., Serial 4159. The rated maximum continuous input of this unit was 460 volt 3 phase, 60 cycle, at 84 amperes and the maximum continuous output of the germanium rectifier was 2000 amperes at 25 volts DC. The sinusoidal power output was obtained by oscillating the auto-transformer input of the rectifier with a four-bar linkage designed by Yang(35) and shown in Figure 31. The linkage was driven by a 2 hp electric motor through a 3/4 hp Vickers hydraulic variable speed reducer and a 10/1 Boston reduction gear. With this linkage an approximately sinusoidal power output was achieved. B. Instrumentation The instrumentation employed is shown schematically in Figure 32 and consists of means for measuring fluid and wall temperatures, fluid mass flow rate, voltage and current.

Figure 30. Heater Exit Fluid Mixing Baffle.

-75Figure 31. Experimental Apparatus Showing Four-Bar Linkage and Power Supply Auto-Transformer.

POWER SUPPLY TO SWITCH PANEL TO SWITCH CURRENT PANEL SHUNT VOLTAGE MIXING fl DIVIDER BAFFLE CALMING II WEIGHT FROM SECTION UI TANK LINE THERMOCOUPLE T T HEATER PROBE WALL THERMOCOUPLES TO SWITCH TO SWITCH PANEL PANEL THERMOCOUPLES VOLTAGE a CUR RENT~ SIGNALS SWITCH PANEL THERMOCOUPLE VISICORDER REFERENCE JUNCTIONS POTENTIOMETER Figure 52. Schematic Diagram of Instrumentation.

-771. Thermocouples Copper-constantan thermocouples were used to measure both fluid and wall temperatures. The wall thermocouples were 30 gage wire spotwelded directly to the outer surface of the test section by a capacitordischarge welding unit. Thermocouples were located on both the upper and lower surfaces of the tube at six inch intervals; a typical installation appearing in Figure 33. The wires were bent in a U shape at the junction with the lead-in end taped to the tube wall after spot-welding the junction. The wires from any one hot junction were oriented one upstream and the other downstream from the junction to minimize the effect.of heat conduction along the wires away from the junction. At the hot junctions a spacing of approximately 1/16 inch was allowed between the copper and constantan wires; the wall material between the wires acting as a neutral metal. The actual temperature measured was the average temperature of the wall between the two wires. When the individual wires were welded to the tube wall, optimum values of capacitor charge and voltage were determined which gave sound welds with little or no weld beads at the junction. This minimized the heat capacity of the junction. With the hot junctions connected directly to the tube surface a DC pickup was observed whenever a potential was applied across the test section for heat generation purposes. This effect was minimized by alligning the copper and constantan wires at the junction as close as possible to an equipotential line on the tube. A calibration of the individual junction DC pickups was made; this appears in Figure 34. This calibration was performed by suddenly applying a voltage across

-78Figure 33. Typical Installation of Wall Thermocouples.

-79+4.0 + 3.0 IU + 2.0 6U I+I.O)- ~~ 5U 1.0 2.0 3.04.05.0O2L93U c3L w 6L a. 2U,4t -2.0 -3.0 4L -4.0 1.0 2.0 3.0 4.0 5.0 6.C0 TEST SECTION POTENTIAL (VOLTS) Figure 34. Thermocouple DC Pickup Calibration.

-80the test section and observing the deflection of the visicorder galvanometer, then suddenly removing the voltage and observing an opposite galvanometer deflection. Since these voltage changes were not steps, extrapolation of the galvanometer trace back to zero time at both application and removal of voltage was necessary to determine the proper galavanometer deflection at zero time. An average value of the deflections observed upon application and removal of voltage was used. In addition, when the system was run under steady state conditions the power leads to the test section were reversed, reversing the polarity of the DC pickup. The difference between the two thermocouple signals observed, divided by two, was the DC pickup. This calibration was within three per cent of the step change calibration for all the wall thermocouples. A Ceramco 22 gage stainless steel shielded, copper-constantan thermocouple was used as a probe at the exit of the line heater where the test section inlet temperature was measured. Both the 30 gage and the 22 gage thermocouples were calibrated at the steam and tin pointso 2o Voltage and Current The voltage applied across the test section was measured by means of a resistance bridge with an output of 2.53 millivolts/applied volt. Current was measured with a shunt resistance having an output of ~00497 millivolt/ampere. 3. Temperature Recording A series of double-pole, double-throw copper knife switches were arranged so that the signals from the thermocouples could either

be displayed on a precision potentiometer (Leeds and Northrup model 8662) or recorded on a recording galvanometer (Minneapolis-Honeywell Visicorder Model 1012). All steady state emfs, voltages and currents were measured with the precision potentiometers while the visicorder was used to record the steady-periodic data. The galvanometers used were Heiland model M-40-120 with a nominal coil resistance of 21.4 ohms, to be used in series with a damping resistance of 120 ohms. With this combination a nominal sensitivity of 1.2 inches deflection/millivolt was possible. In order to improve the galvanometer sensitivity the damping resistors were removed and the galvanometers were driven undamped, except for the inherent resistance of the thermocouple circuit, which gave a sensitivity of from 1.48 to 1.82 inches/mv, depending upon the thermocouple resistance. With this difference in the external load resistances there theoretically was a slight phase shift between galvanometer traces. However, the frequencies of the signals recorded in the experiments were all less than 2% of the natural frequency of the galvanometer (40 cps) so that the phase lag was very small and was not discernible on the visicorder trace. A phase lag of the same order of magnitude existed between the galvanometer recording the test section voltage and the galvanometers recording thermocouple emfs. All of the thermocouples on the upper centerline used the same reference junction; one thermocouple at a time being connected to -the junction by a knife switch. The same type system was used for the thermocouples on the lower side of the tube. The heater exit thermocouple probe had its own reference junction as did another entrance

-82thermocouple. This last thermocouple was spot-welded to the test section four inches upstream from the first power electrode. The purpose of this thermocouple was to act as a backup for the heat exit thermocouple probe and was also used to determine if there were any appreciable heat conduction upstream through the tube, away from the heated section. 4. Mass Flow Rate Fluid mass flow rate was measured by means of a weigh tank at the exit of the test section. An average of five flow measurements was taken at the start and at the completion of every test run. C. Test Procedure The experimental program was divided into two major parts: steady state and steady-periodic state. For both the steady and steadyperiodic states the Reynolds number was varied from 1000 to 2200 with mean heat generation rate being varied from 100,000 to 450,000 BTU/hr-ft3. For the steady-periodic state the period of the heat generation was varied from 5 to 20 seconds with a range of E from 0.1 up to 0.8. For any one test run the water was turned on and the inlet temperature adjusted by the line heater. Mass flow rate was roughly adjusted to give the desired Reynolds number based on entrance temperature. Power to the test section was adjusted to the desired level. For the steady state this was simply accomplished by remotely adjusting the input to the power supply autotransformer. Steady-periodic power level was set by manually adjusting the autotransformer shaft to the mean power level, adjusting the oscillation linkage to give the desired C and setting the

-83output level of the hydraulic speed reducer to produce the heat generation frequency for that run. When the system had reached either the steady or steady-periodic state (about five minutes) the mass flow rate was again adjusted to give the desired Reynolds number based on the mean steady mixed-mean temperature, which is given by ( \)AV 7- ~ 3.4/13 WATr7s After another waiting period, temperature measurements were made. D. Data Reduction With the thermocouple signals having the DC pickup superimposed, an involved data reduction procedure was necessary. For each of the six stations along the axis of the tube a Visicorder trace was obtained which contained voltage divider and current shunt emfs, upper and lower thermocouple emfs, and the two entrance thermocouple emfs; a typical visicorder trace appears in Figure 35. At each one second time line on the trace the following reduction procedure was followed. Upper and lower wall thermocouple galvanometer deflections were measured and converted to my. The voltage divider and current shunt galvanometer deflections were also measured and converted to test section volts and amps. After computing the power delivered to the test section, the DC pickup calibrations were multiplied by the applied voltage and the thermocouple emf's corrected for DC pickup. The resulting numbers, the thermoelectric emfs, were corrected to standard and converted to ~F. The same procedure was followed for the entrance thermocouple emfs

Figure 55. Typical Visicorder Trace, during Steady-Periodic Operation. A. Applied Test Section Voltage B. Test Section Current C. Upper Wall TIC D. Lower Wall T/C E. Test Section Entrance T/C F. Heat Exit T/C

-85with the exception that these signals had no DC pickup. Subtracting the entrance temperature from the wall temperatures gave the ATs for the upper and lower wall. The factor was obtained from the kf equation k,,. pR L (Ro'- R) The dimensionless temperatures at the upper and lower wall were then obtained by dividing the ATs by this factor. An arithmetic average of the upper and lower R*(2) was taken to give for that w Ri ine. This procedure was AVG time line. This procedure was repeated at each one second time line until a full cycle of applied voltage was covered. The applied voltage was then plotted along with the dimensionless temperature response. A maximum and minimum temperature were determined from this plot, from which were determined (cA% Aoe = rted o Vrain- (s)M 4 [ ( *4tAX Ow ( *vMI The phase lag of the temperature response was measured in seconds and converted to radians.

-86Eo Results of Experiments During the preliminary test runs it was discovered that free convection effects were present to such an extent that they completely dominated the forced convection except under the limiting conditions of large laminar Reynolds numbers and small energy generation rates. The remainder of the tests were run under these restricted conditions and a small difference still existed between the temperatures measured on the upper and lower surfaces of the tube. With small energy generation rates the temperature differences existing in the system were small and the temperature dependency of the physical properties could not appreciably affect the data. Figure 36 presents the steady data representing the arithmetic average of the measured steady temperature at the upper and lower surfaces of the tube at r - R- as a function of axial distance. The data Ri correlate well with the numerical solution. At the entrance the data show a tendency to diverge from the numerical solution. This is due to the fact that the entrance condition applied in the computer solution in the wall is idealized. It is impossible to have heat generation in the wall at x = 0 and, at the same time, to have zero temperature for all r there. In fact, the steady wall temperature difference between r = -- and r = 1o0 must be invariant with x. Thus at the entrance Ri the temperature at the outer perimeter of the tube is greater than zero but falls rapidly to zero upstream from the entrance, as shown schematically in Figure 37. The temperature at any radius in the wall will decay exponentially in the -x direction, upstream from the entrance, due to conduction.

-873.8 3.6 3.4 3.2 3.0 - 2.8 2.6 2.4 - 2.2 O O 2.0.8 8:, -8 1.6 1.4 1.2 1.0 NUMERICAL SOLUTION ~~0.8 0 ~O EXPERIMENTAL 0.6' z.O43 8:0.126 A 0.4 0.2 0.002.004.006 008.010.012 RX( RePr) Figure 36. Steady Wall Temperature at r = Ro/Ri versus Axial Distance.

-88IDEALIZED OUTER WALL TEMPERATURE *sw CTUAL OUTER WALL TEMPERATURE 0 O * ( Re Pr)-' Ri tsw ( r= ii q *sw (X T' x Ri Ro Figure 37. Schematic Diagram of the Actual Wall Temperature Distribution at the Entrance of the Tube.

-89In Figure 36 the majority of the data lie above the numerical solution. Although some of this deviation is due to normal scatter, a small increase in temperature is due to a film layer deposited by the water flowing in the circuit. This film was removed periodically but reformed almost immediately after flow was reinstated. Figure 38 presents the measured steady-periodic temperature at the outside of the tube as a function of time, compared with the numerical solution corresponding to the experimental conditions. The data fall, in general, slightly above the numerical solution and show a slight phase shift. The test section tubing being relatively thin walled (8 = 0.126), the computer solution used for both Figures 36 and 38 was the lumped-wall program. The analysis in Appendix D shows that the steady wall temperature difference for the conditions in Figure 36 is 3.32 x 10-4 or approximately two per cent of the predicted outer wall temperature at x =.001. For the steady periodic data in Figure 38 the analysis shows a comparable temperature difference between the inside and outside of the wall with a phase shift of less than one half a degree between the two temperatures.

180 0 q,, |\NUMERICAL 140_ \ O EXPERIMENTAL [POWER] (WATTS) 100 60 0~~~ 0 20 ~ 0.793 2.8 0 r' 216 & a 0.126 4, (\ R O1.16 a _....I I 0 2 4 6 8 10 12 t ( SECONDS) Figure 38. Steady-Periodic Wall Temperature at r R/Ri versus Time.

CHAPTER V CONCLUSIONS 1. In formulating the basic differential equations governing heat transfer to laminar flow in tubes, the assumptions which are made may not be realistic in that the regions in which they are valid may be quite narrow. In particular, the effects of free convection can become significant and, in the limit of small laminar Reynolds numbers and large energy generation rates in the tube wall, free convection can be of such a magnitude that any physical model neglecting it would be erroneous. 2. The periodic Nusselt number which is observed when the energy generated in the wall is represented by the sum of a steady component and a sinusoidally oscillating component is itself the sum of a steady and a periodic component. The time average of this periodic Nusselt number is equal to the Nusselt number representing the response of the system to the steady component of energy generation in the wall. 3. Any extension of the classical Graetz problem which includes transient or periodic conditions should include the influence of the wall heat capacity. The wall heat capacity has a marked effect on the steady-periodic response of the fluid temperature adjacent to the wall and it is felt that the same would be true in the transient state. 4. The effect of generation frequency on the temperature response is an increase in phase angle at all locations in the wall and fluid as frequency increases. The maximum phase lag in the wall is - radians 2 while phase lag in the fluid increases continuously as frequency increases. -91

-92The limit of the periodic amplitude of the temperature response in both wall and fluid is zero as frequency approaches infinity. In the other extreme, the limit of the periodic response of the system as generation frequency approaches zero is the steady response of the system to the instantaneous value of energy generation. 5. The problem simplification gained by assuming a radially lumped wall is justified when the wall is thin and when the ratio of fluid thermal conductivity to wall thermal conductivity is small. 6o A physical model which assumes a zero entrance condition in the wall when the energy generation is spatially uniform is in error. A more accurate model would perhaps include a wall temperature distribution at the entrance obtained from a simplified but independent model.

APPENDIMCES APPENDIX A DERIVATION OF THE DIFFERENTIAL EQUATIONS The system under consideration, shown in Figure 1, consists of a constant diameter tube through which a fluid moves in laminar flow. Heat is generated in the tube wall at a uniform volumetric rate which varies sinusoidally with time according to the expression 9-= 9 (0 + n ~nt) (A-1) Heat is transferred from the tube to the fluid while the outer perimeter of the tube is adiabatic. The first law of thermodynamics, within the framework of the assumptions made in Chapter II, applied to a control volume in the fluid, Figure 39, yields pys a 4- CzLJ3a~S~ 3T~OfC (r)as (A-2) The conservation of momentum applied to the same control volume gives the radial dependency of the fluid velocity u1 (r*) - (/- R; 12 ) (A-3) where Umax dp Ri the centerline velocity. Further, a mass baldx* 4t ance across the entire fluid cross section yields UAx = U-93-93

-94rq +dr* q r r U(r)'.,, /////' dr" rr q q x A01- g dx" --— Ri Row \~~~~~~ \ -x* 4- dx —-- l --................. Figure 39. Control volume in the Field. sqr+dr' _L I] d~~ry —.-. x - d *.,> R< _~~~~ re V'.////////////; Figure 40. Systemn in the Wall.

-95where Um is the average velocity. Equation (A-2) now becomes 2 - )A(A-4) The first law of thermodynamics applied to a system in the wall, Figure 40, gives Nt ur dr a r) *'L Cs (l - + Sitn W) (A-5) Equations (A-4) and (A-5) are coupled at the fluid-wall interface, r* = Ri, with the equality of the fluid and wall temperatures there. In addition, the solution of Equation (A-4) must satisfy the condition 9r ('c, t,o)=o for reasons of symmetry, and the solution of Equation (A-5) must satisfy the boundary condition due to the adiabatic condition imposed there. Equations (A-4) and (A-5) are also subject to the entrance condition T(o,t,r*) =To Equations (A-4) and (A-5) can now be non-dimensionalized with reference to the following set of dimensionless variables and parameters:

-96go Rio Re 2 L- R= ttc 2 ~to _ ~~- A r,Q, Equation (A-4) becomes a(L r~_,3 a(LE _ al~Lf ~i- r3L ~ O C (A-6) se' % ar~ r ar and Equation (A-5) becomes (y &r )2Of(/+ 6 nTO) XL~ (A7) which must satisfy the conditions r (%Q) 0 =r('L B. R ) -

APPENDIX B DERIVATION OF THE FINITE DIFFERENCE EQUATIONS 1. Coordinate Notation The fluid and wall are divided up into a network both radially and axially; the temperature at the center of mass of each concentric element representing the average temperature of the material in the element. This division of the two regions is shown in Figure 41. In referring to the location of an element in the fluid, the subscript (I,J) is used; I representing the Ith element measured radially from the centerline of the tube and J representing the Jth element measured in the positive axial direction from the entrance of the tube. A point represented by (O,J) is on the centerline of the tube while a point (I,O) is at the entrance. In the wall region a point (K,J) represents the Kth element measured radially from the inner wall of the tube and the Jth element from the entrance. A point in the wall represented by (O,J) is a point at the fluid-wall interface while a point (K,O) is at the entrance. The spacing between the elements depends upon the number of elements employed in the network - sometimes referred to as the "grid size". If there are N elements between the centerline and the inside wall of the tube (not including the element on the centerline) the dimensionless Ar is 1 since the dimensionless inner radius of the tube N Ro wall is 1.0. If the dimensionless wall thickness (Ri - 1) is 5 and there are M elements radially in the wall (not including the element at the fluid-wall interface) then Ar in the wall is 5/M. The radius of -97

o( KJ) o(K,J+I )I I I'I.. I - I -I-_ - Figure 41. Wall and Fluid Grid Network. I I I I 0r I- -- ~ - -- -- ~ - I 0 0 1 I- (f, J) I (I, + I --

-99an Ith point in the fluid is I/N, (I = 0, 1, 2...N), and the radius of a Kth point in the wall is 1 + K- (K = 0, 1, 2...M). In both the fluid and the wall regions, for a tube of length L divided into P elements (not including the one at the entrance), the dimensionless Ax is L/P and the dimensionless axial distance from the entrance to the Jth element is JL. Thus, the coordinates of the point (I,J) in the fluid are (I/N, oJL) and for the point (K,J) in the wall the coordinates are (i + K, JL). M P 2. Finite Difference Derivatives When attempting to write a differential equation in finite difference form an error is introduced whose order of magnitude may be estimated with the aid of Taylor's series. The Taylor series of a function f(x,y) at a point (x + Ax, y) about the point (x,y) is fI(+ ) ()) (%) f+ + ( as (B-l) and the function at (x - Ax, y) is expressed by JI 4) + _ _. (B-2) Subtracting (B-2) from (B-l) yields SI(~C~ ts - t,) 2 A)C +23... so that (&f) - As >L - %,if) (B-3)

-100where the order of magnitude of the neglected terms is Asx2 when Nx is small. In a similar fashion, when (B-i) is added to (B-2) cf3 = f(%LL, ) L (- %, ) - (%,) + L (AL)) (B-4) Again the order of magnitude of the neglected terms is nx2 when Ax is small. In some instances (when the value of the function is unknown at (x + Ax, y)) the derivative of the function at x can be expressed in terms of the values of the function at (x,y) and (x - Ax, y) by rearranging Equation (B-2) __,,.~, AM(B-5) where the order of magnitude of the neglected terms is Ax. Equation (B-3) is thus a far better representation of the derivative than is (B-5) when Ax is finite and Equation (B-2) is used whenever possible, 3. General Steady Finite Difference Equations In Chapter II the function V was divided into steady and periodic parts; 4s and *p. The differential equations for the steady functions Asf and tsw were j(-r ) ~a~i'a4~ X_ 4- (Fluid) a~c r~' r ar (B-6) a 4- / DISL = o (Wall) a I-" r

-101 In the finite difference representation of the derivatives at a general interior point (I,J) in the fluid or (K,J) in the wall the following notation is used Fluid Wall ICISS (,J ) S 5(I) Y )S (K, Jo- /(1) =, (:-J J') v S')', (:- (K-I) VS (Iz, z-/): _ U) These one dimensional arrays are more versatile in the computer program than would be a two dimensional array. A marching solution is used in the axial direction in which a value of x, i.e., J, is set and I is varied from O to N in the fluid and K from O to M in the wall. At each successive value of I or K the difference equations for that point are solved in terms of the current values of the variables at the surrounding points which make up the independent variables in the equations being solved. If the stability requirements discussed in Chapter II are met, all of the variables at this value of x will converge to their true values if iteration is continued long enough. The representation of the derivative 6sf in Equation (B-6) ax must be the one given by Equation (B-5) since Vsf(I, J+l) is still unknown when computing the value of Isf(I,J). Thus _ s(I_) - I () s (I)-U(I) (B-7) L3L n~1 L/

-102Using (B-4) to represent the second partial derivative of *sf with respect to r gives ar-e ~(sI f5(-1)- (sVi) (B-8) ar~ The last term in the steady fluid equation is, using Equation (B-3) to represent the derivative, __ r /N ___~ _(B-9) r ar _ /(;/ S )) Substituting Equations (B-7, -8 and -9) into the differential Equation (B-6) and solving for S(I) gives Ns)( N(I i) s(/) 4- Ns (/( ) - ) - () u) 2N'+, - (_B-10) This equation is for a general interior point and is not valid for either I = 0 or I = N. Inspection of Figure 41 shows that these elements are not full elements and so special equations are necessary for them. The derivation of these equations requires a return to the basic problem. This appears later in this Appendix, The equation for a general interior point (K,J) in the wall, using Equation (B-8) and (B-9) in the wall differential Equation (B-6), is V ( A)i t Vk (-) -2 V(g) i (V(K,,)- V(K —I) - o

-103Rearranging and solving for V(K) gives -( M ( 12(B-11) Here also the equation is valid for a general interior point in the wall but is not valid for K = 0 or K = M. The equations for V(o) and V(M) are derived later in this Appendix. 4. General Periodic Finite Difference Equations The periodic functions *pf and. pw defined in Chapter II must satisfy the equations a e is art r ar (B-12) a_ _ _ _a. 4- F-r ) + c SLnre Ka I r ar It was assumed that the periodic functions would be of the form ",pf = A r, ) SinLe + B(x,r) CosIre (B-13) and /,'/,u ~ C(%,r) Sin re D(x,r) Cosre (B-14)

Substitution of (B-13) into (B-12) yields AICosre - BrsPnre + [/- -A Sinre + Cost = (r 9A r ar I B> osr Gar" r- ~ os r r This equation is valid for all @ ( >> 0)o One can combine all Sin f9 and Cos'G terms and equate them to O. Thus Al +(-r ) a a_ B, B - B_) ar+ r ar -TlrA LA + /A a rL r ar These two equations can now be written in finite difference form for a general interior point (I,J) in the fluid using the following notation A(I,r) = A(r) B(,r) = s() A(z-IJ) = A (z+I) S B(I+I,,)- B(z-l) A(IJ-I) = (I)(IJI) F() Thus.., tX — B(r) - hl',-.,- x. -..,14 4, B (I-, r -PA +P /-(~N)

-105Rearranging and solving for A(I) A' N ( r); B(= Z A ) (4- P)) B(I-) -j ~~ ~~~~~2NZ~ ~~~ (z i- ( - (B-156) (r Lr ( Lr N By similar substitution of Equation (B-) into the wall differential Equation B-12)for B(I) is gives JTC Cos DO -JD Sair O c / ( e nC*' s nZ'r +.,F + r.. + Separating the Sin PG and the Cos PG terms as before,C =A (i. / __p -rc = t O r r 8r' rD=iffe rence fo rm the equation for is In finite difference form the equation for C is JTC(K)= _ jD(f/l)4+ D(K-4)-2D(k) / (Dfr+,) - A (Y) (I4 /< Z J

-106Upon rearrangement C I~) = t ( + 2 D{~+e) + J~ (M ___ N (B-17) - 2/. M~ D(t) Similarly A-6 M / (B-18) t z Mz c(K)- ~ r I& I The periodic functions at r = 0, 1, and Ri require special attention Ri similar to the steady functions at these points. 5. Limiting Elements and. Their Equations ao Centerline Element The boundary conditions imposed at the centerline were a s,, (,V, o) a (a,-( )', 0~) - o ar ar The last term in the fluid Equations (B-6) and (B-12) is 1 r ar which is undefined, at r = 0 since Lirn / a-/ o r — or r o

-107It follows from l'Hospital's rule that Li rn r;5r C- r -~ Thus, at the centerline of the tube, the steady and periodic Equations (B-6) and (B-12) respectively are (4IS =, a _s) (B-19) a?C 3 r +P 2 (B-20) Substituting (B-4) and (B-5) into (B-l9) and solving for S(o) yields ~5(o) = +j\(I~ () / LI) (B-21) 4N~, 914_ Substituting (B-13) into (B-20) and equating Sin Ir and Cos PG terms to 0 as before gives the two equations - BT + A = 52 aL In finite difference form these become A(o) = 4rN B(l) - 4r ) 8(o) P F(o) (B-22) ) r Lr z~ ~ ~ ~ ~ ~ ~

-108The same equations for S(O), A(O) and B(O) occur when one returns to the basic energy equation and derives the steady and periodic functions from basic principles. b. Fluid-Wall Interface Element The point (N,J) in the fluid and (O,J) in the wall are identical; the point representing the last (radially) 1/2 element in the fluid and the first 1/2 element in the wall, as shown in Figure 42. With the element being composed of both fluid and wall material, it is necessary to return to the energy equation to derive the element equations. The energy equation for this composite element is (uAurC + (OAf ) 7[ - Ain - ut Aof 9 Ur where Au 2Z7rR& 6, A S - 2 kk r*-r\ 40L, 2 - wP", 90L~ R In terms of the previously defined dimensionless variables I, r and @ and the parameters r1,, and A this becomes (, a = MRS _ a) \Mf/ e6 R. I (B-24)

-o109qlN - -,q It (IN WALL) (I = N J ) (K=O' Rf Rw \~~~~~~~~ d- r Fiture 42. Fluid-Wall Interface Grid Element.

-110Again s = ~s + *p, which yields the two equations RS 4- Sur R (B-25) \r RL R rr #? R. Ir =I R4 ( / + A 4 d~k- Si(B-26) e P 3r/r= Ri r, IZ From the defined grid spacing in the wall and fluid d,.. R,_ / +.R, 2 M Re 2, 4,_ R7 _ /- / R, 2N R 2ZN Employing these dimensionless distances and radii and the finite difference derivative notation, Equation (B-25) becomes _/ _( _ (/ I __ _ S(N\)- S (N-I) = O I V / IN/ 2M where V(O) = S(N). Solving for V(O): V(o) #= (7 4)V(i) 4 (NA)s(/) 2/ (B-27) o (3 n (2) (- 2) Substituting Equations (B13) and (B-14) into (B-26), separating the coefficients of sine and cosine and employing the finite difference

derivatives yields - /r 2' 2/") D(o) =.A. -C (I) -) C (o) _ A (N) - A (N-1) ) + V/N Z z, and 2M 2 B (N)- B(Nl _ / Solving for D(O) and C(O) while noting that C(O) = A(N) and D(O) = B(N) C(o) = __ (>j ~ )DO) (- M + +)4N-i)Do (B-28) E o N- _-__N-_ _ -(NU-)A (N- I) - } h(B-29) c. Element at the Outside of the Wall Referring to Figure 43, the point (M,J) at the outer surface of the wall represents an element which is one half of an element in

-112q"(I1N WALL ) -(MAJ qOUT A Rw Fiue 3 Gi leet tth usieofte aI! I NSULATION dx Figure 43.. Grid Element at the Outside of the Wall.

-113thickness and whose outer surface is adiabatic. The energy equations is d T7 rAI C' a AorCt a 9-Au - where Aw: 27Ro.6,: U k Tu A,: 27kRur In terms of the dimensionless variables and parameters this is R4 -6+ (i (B-30) R. RX as RL ( r R=- S Dividing?w into its steady and periodic components and further dividing the periodic component into sine and cosine terms yields the three equations FRo C - A Rur where RR =/ Rc = /: 2 E

-114If we now write these three equations in difference form and solve for V(M), D(M) and C(M), respectively V(M) = V(M-I) (/+ t) z4A M (B-31) 2/i (7 )(B-32) -a( M M-X) D (M) + A ( + M- D (M-1) C(M)= (B-33) -(", /.a) 60 Lumped Wall Equations For the case where the wall is very thin, a radially "lumped" wall approach is necessary. Here it is assumed there are little or no variations in wall temperature in the radial direction. Referring to Figure 44, the temperature at the fluid-wall interface represents the temperature of the last half element in the fluid and the entire wall. The energy equation is (I - A,C,+ /sA, C) -_ Aout go Ar9 where A- -'7R,-(4-P;R Ao.t = 2wRr Af -- 27rei &e

-115INSUL ATION ~WALL q v FLUID (NJ) qOUT Ri Rf R f dx x f~ Figure 44. Fluid-Wall Interface Grid Element for a Radially Lumped Wall.

-116Employing the dimensionless variables and parameters previously defined, this becomes R / 1 2 ae 3 z - -2 R r I 8 $S V {6 V R, / eJ;59 R t r re (B-34) The dimensionless thickness and radii are 4;~ = I R. I' R2 / = $~ Re 2 N R N R. Equation (B-34) now becomes ( 26 =-). r= +2)(Is r) Once again, s = ~s + ~p such that (2- ()( ) 4 (8 26 ) =0 ~ (B-35) ( 26 4) ( 2- (e ) (6 2&SLnPre (B-36) Equation (B-35) in difference form is s (lN) = S(N- ). \ 2+/I (B-37)

-117and the periodic function, divided again into sin and cosine equations, gives the two difference equations. A(N) - (2- )B(Q) c (2N-) B(N-I) (B-38),(,.. / +-. "-t~ 2.~ B

APPENDIX C COMPUTER PROGRAMS The finite difference equations listed in Table I, valid when r is large, are solved by iteration in the first program listed in this Appendix. Those equations listed in Table III, valid for small r, are solved in the second program. In both programs the nomenclature used was kept as close as possible to the nomenclature of the text. Those special functions and redefinitions which were necessary are listed in the computer nomenclature. In both programs the following sequence of operations is followed: read in data, test for stability, compute constants, solve difference equations, compute phase angles and amplitude ratios, compute Nusselt numbers, print results, increment x, solve difference equations etc. This sequence is shown in Figure 45, the flow diagram for both programs. -118

-119EAD DATA PRINT DATA PRIN1T DATA IS DATA REJECT STABLE F STATEMENT COMPUTE CONSTANTS INITIALIZE VARIABLES SOLVE DI FFERENCE EQUATIONS HAVE ALL IARIABLES CONVERGED/- I INCREMENT PRINT DATAI t T x REJECT T S CONVERGENCE SATEMENT TOO SLOW COMPUTE IS L-P AND AR T COMPUTE PRINT Nu|, Nup, -- RESULTS Figure 45. Computer Program Flow Diagram.

-120PROGRAM A - LARGE r DIMENSION A(25),B(25),C(25),D(25)E(25) F(25) G(25)H(25) 1Q(25),R(25),S(25),U(25),V(25),W(25)tY(25),PHIF(25),PHIW(25) 2PERCOF(25),PERCOW(25), EPSW(25),OPPF(25),OPPW(25), 3ADJF(25),ADJW(25),AMPRAF(25),AMPRAW(25) INTEGER I,J,KM,N,P,COUNT PRINT COMtIENT$2$ PRINT FORMAT TITLE START READ DATA PRINT COMMENT $1NEW DATA$ PRINT RESULTS EPSQ,BETA,LAMGAM,DELTA#LtP DELT = 10.OE-37 ALPHA =.001 PI = 3.14159 THROUGH LOOPTFOR N=20,-1,N.L. 10 WHENEVER GAM *G. 8.*(N.P.2.)+P/L PRINT COMMENT$OTHE NUMBER OF NODAL POINTS IN THE FLUID IS$ PRINT RESULTS N TRANSFER TO ENTRYA OTHERWISE CONTINUE LOOPT END OF CONDITIONAL PRINT COMMENT$ODIVERGENCE OCCURS IN THE FLUID EQUATIONS WHEN 1N IS REDUCED TO 10. DATA IS REJECTED $ TRANSFER TO START ENTRYA THROUGH LOOPUFOR M=109-1,M.L. 2. WHENEVER GAM.G.4.*LAM*(((M/DELTA).P.2.)+(M/(2.*DELTA))) PRINT COMMENT $OTHE NUMBER OF NODAL POINTS IN THE WALL IS$ PRINT RESULTS M TRANSFER TO ENTRYB OTHERWISE CONTINUE LOOPU END OF CONDITIONAL PRINT COMMENT$ODIVERGENCE OCCURS IN THE WALL EQUATIONS WHEN M 1IS REDUCED TO 2. DATA IS REJECTED $ TRANSFER TO START ENTRYB CA = N.P.2. CB = P/L CC = CA/GAM CD = CB/GAM CE = M/DELTA CF = CE.P.2. CG = CE/2. CH = 2.*CF CI = BETA/LAM CJ = LAM*C(F/GAM CK = LAM*CG/GAM CL = 2.*CJ CM = BETA*EPSQ/GAM CN = 4.*CA CO = CB + CN CP = CN/GAM CQ = CO/GAM CR = CB/GAM CS = (CE + 1./2.)/CI CT = N-1./2. CU = l./(2.*CE) CV = CU*EPSQ CW = (1./(2.*BETA*CE)) CX = (CW+1./(2.*N))*GAM

-12 1CY = CE+M-1./2o -----— i —------ -- - - - - - - - - - - - - - - - - - - - ---- -- ------ DA = CZ*BETA/LAM DC = CZ*BETA*EPSQ DD —=...C-Z 4G'AM DE - 1./N DF = 1/(N.P.3.) PRINT COMMENT $OTHE CONSTANTS GENERATED BY THIS DATA ARES........PRI NT'-R E!SUZT..'~CC'C-q ~T3E C'~-CF~'CG-'"JC_ C K, C L",'t M, C N, C 0, 1CP, CQ,CR, CS t, CT, CU, CV, CW,CX, CY.CZ,iDA, DB, DC, DO DE DF THROUGH LOOPAFOR =0','I.'G.N. A(I) = 0. -' —BTIQ~.... —-- -—,. ——. ——... S(I) = 0. -.........El(I) -"-; --................................... F(I) = O. LOOPA U( I)' = 0 THROUGH LOOPVFOR I=0,1,IoG.M.............. C (I.)...= 0 ~.......................... D(I) = 0. LOOPV V(I) = 0....... THROUGH LOOPB,FOR J=1,1J.G.PP THROUGH L'O-PCF 0-T;11C4 —R —- — I ——.;-;.........G(I) = 0O ----— )-0. ------------------------------ - FF~~~~i-F-y- — C. --- ----- LOOPC W(I) O, THROUGH LO'GP'D~'FOR- OR; o...-.......... Q(I) = 0......................................- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - R(I) = 0. LOOPD Y(I) 0....................... ~6ON —----- ------ --- -- --- --- --- --- -------------- -- COUNT = 0 ENTRYC D(M) = (DB/DD)*C(M)- (DB/DD)*C(M-1)-DC/DD... C ( M')'~ ~- —'... (> bi B'-/'-D r~s 1)"~ D~I'~F`~-3'TD-/lcIc~T) M-2 -'-I 1 -)........ V(M) = V(M-1)+DA/CY THROUGH LOOPE'FOR I=M-1,-1,I.L.1 D(I) = -(CJ+CK/(1i.+I/CE) )*C( I+1)-(.CJ-CK/(1i+I/CE))*C(I-i)......... IzFCE:~-C.-Ci.T —- — C —---- C(I) = (CJ+CK/(1.+I/CEI))*D(I+1)+(CJ-CK/(1+I/CE))*D(I-i) 1-CL*D (I) FAC = 1./(4.*C CE+I)) -LOO P0E VNTV-'-=' ( (T.T2 +F AC *iV(I+1) +T ((1 * /2.-]fA) *V/(-1). +CI /CH —-- D(O) = -(CS/CX)*C(1)+ ((CS+CT)/CX)*C(O)-CCT/CX)*A(N-1)-(CV/CX) --— z —------— l- ) V(O) = (CS/(CS+CT))*V(1)+('CT/(CS+CT))*S(N-1)+CU/ (.CS+CT). -AIN ) = CC ) 8(N) = D(0) S(N) = V(O) THROUGH LOOPFtFOR I=N-l,-1I.L.1 A ( I) -* ( -*+1. / (2* I ))B ( I+ 1)i+CC* (!.- 1. /.(2.*.I.) ) *B ( I-1)-(CD 1( 1.-( I.P.2. )/CA)+2.'*CC)*B ( I) +CD*( 1,-( I.P2, )/CA)*F( I) B TTC ) —:'T:-CCI- +T./(2.~'I) )1*A ( I1 +*l-.T-CC' I,-l,/'(2.*I ) i)*,A (I-I) + (CD i*(li-(I.p.2 ~)/CA )2+2 *CC*A ( I)-CD* (,-( IPc2)/CA)*E (I)..L0~ — PF. —---— (TI- CA-K* I -1 1 - 2+. 7T./ C('*I). /CB* C 1.- C P2.2~ ) /CA )+2,*CA) )*SA (I 1+(CA*( 1.-1. /(2,*I1))/ (CB*(l1,-(C'2. )/CA)'+2.*CA) )*S(I-1)'2-f —-C~-i- -I-; —-Pg-~-FT-CA') / (CB*( 1,~~TP-( I S-F~Sr 2 ~)/CA) +2,*CA)')*CU( A l) A(O) = 4.*CC*B(1)-(4.*CC+CD)*B(O)+CD*F(O) -- (0) = -4.*CC*A(i )+z(4.*CC+CD)*A(O) -CD*E-(O) S(O) = (4.*CA*S( 1)+CB*U(O )-/(/ +.*CA +CB)

-122THROUGH LOOPGFOR I=0,1,I.GeM WHENEVER.ABS.(C(I)).L. DELT.OR..ABS.(D(I).L.'DELT.OR..ABS. 1(V(I))..L.DELT TR ANSFER TO - ENTRYD OTHERWISE ------------—' —-G —— EWH VER —,-;-TBT(C' T-';QtC-T'7 C ( I ) ).L*ALPHA.AND.ABS. (D ( I)1R(I))/D(I)).L.ALPHA.AND..ABS. ((V(I)-Y(I))/V(I)),.L.ALPHA --------------— CTINUE —------ ---- OTHERWISE WHENEVER COUNT *G. 500 -...............PRINT COMMENT $ ~*:*******************************************$ PRINT COMMENT$ CONVERGENCE IS TOO SLOW. DATA IS REJECTED$ ----------------- rNrTP-cOT-'fET$-CO $ -SE-OVc F- ONVERGEN-CE IN —- WALE:-Z'EQUAT'I ONS * $ PRINT RESULTS S(0)..S(.N),W(~0)...W(N),A(O)-.*A(N),G(O)...G(N) Ii BTO —)''FTF.H N)-, V O ).~.V(M), Y(O).. Y(M), C( 0) C (M), 2Q(.0)..Q(M),D(0)..D(M),R(O) R(M) --- ~.... PRIrNT.-..COMMENT' —$-}- ~-*-* *~* -***' * * * *'** * *'*'* * * ** ***** **** ** ***** ***$ TRANSFER TO START ----------------— THERWTSE —-------- ------------------------ CONTINUE END OF CONDITIONAL THROUGH LOOPHFOR K=O,1,K.G.M.....-''(K.) = CCK'.-l-. —---- R(K.) = D(K) — LOOF — Y VK —---------------------------—. —-------------- THROUGH LOOPI,FOR K=O,1,K.GN.-...... G'(K) ='A(KV'. H(K) = B(K) — C'OO'PI..- -.W-(' K —) —''.- S-KT..-. — TRANSFER TO ENTRYC ----—. —---- END OF- CONDITIONAL ----- ---- ENTRYD END OF CONDITIONAL LOOPG CONTINUE THROUGH LOOPJ,FOR I=O,1,I.G.N ------------- -LW-F E-NEVER.ABS'A-'-).. —ELT.O-,-'T-RA' B S(B )" -.'.D ELT'OR.ABS 1(S(I)) L.DELT - -TRA-NS-FER. — TO- -ENTRYE --------... —---—....- -- OTHERWISE WHENEVER.ABS,((A(I)- G(I))/A(II))LALPHAAND,'ABS.((B(I)1HI)) /B(I)).L.ALPHA.AND*,ABS, ((S(I)-W(I))/S(I))*L.ALPHA...................CONTINUE-............................... —------ OTHERWISE ----------------— C — ---------- WHENEVER COUNT.G. 500 PRI-NTCOMENT — E -T-sF'-,-'-.-, *".*",-* ~,',' * * * *** ** $ PRINT COMMENT$ CONVERGENCE IS TOO SLOW. DATA IS REJECTED$ ----------------— PRINT —COMENT —-OW — CONVER-G-ENCE — I'N — FL U I D EQUATI ONS * $ PRINT RESULTS S(O)...S(N),W(O).o.W(N),A(O).*.A(N),G(O)*oG(' N) 1iO —-----------— N —'FF-(;-;H-N),V(O)...V(M), Y(O),,.YiM), C()..C (M), 2Q(O)...Q(M),D(O)...D(M),R(O)*..R(M).PRINT COMMENt $- " —s*- *-*********4*************************-* TRANSFER TO START ----------------— OTHERWISE —--------- ----------- ------ CONTINUE -N —W —----------— ElMO-OF-CONDTTTCAL —-- —.-. -- - -.....- -' THROUGH LOOPK,FOR. KO,1,K.GsM Q(K) C D(K) R(K) D( (........

LOOPK Y(K) = V(K)...'-...-.THRRU "-CE-E'iFO` FI,'iiG'.N G(K) = A(Ki)...H'R(I'r-"B K. LOOPL W(K) = S(K) TRANSFER TO ENTRYCEND OF CONDITIONAL ENTRYE END OF CONDITIONAL LOOPJ CONTINUE THROUGH TOUPM —FR....I=O,TiTI-;G*N E(I) = A(I) F(I) = B(I') LOOPM U(I) = S(I) X = J/CBPRINT COMMENT $0$ PRINT FORMAT XUTX - PRINT FORMAT ITTAB,.COUNT THROUGH LOOPN,FOR I=O,l1,I.G'N NORMF = 1/(.ABS.(B(I))) OPPF(I) = -B(I)*NO'RM ADJF(I) = A(I)*NORMF PHIF(I) = ATNi.(OPPF('I),ADJF(!)) LOOPN PERCOF(I) = (SQRT.((OPPF(I).P.2o)+(ADJF(I).Po2.)))/NORMF THROUGH LOOPO,FOR K=O,1,KoGeM NORMW = 1i/(.ABS~(D(K))) OPPW(K) = -D(K)*NORMW ADJW(K) = C(K)*NORMW PHIW(K) = ATNI.(OPPW(K),ADJW(K)) LOOPO PERCOW(K) = (SQRT.((OPPW(K)oP.2o)+(ADJW(K)oP.2e)))/NORMW THROUGH LOOPQFOR I=O,1,I.G.N LOOPQ AMPRAF(I) = PERCOF(I)/(S(I)*EPSQ) THROUGH LOOPR,FOR I=O,1,I.G.M LOOPR AMPRAW(I) = PERCOW(I)/(V(I)*EPSQ) SUMS = S(O)/(2.*CA) SUMA = A(O)/(2.*CA.) SUMB = B(O)/(2.*CA) THROUGH LOOPP,FOR I=l,l,I.G.N-l RADFAC = (I*DE-(IoPo3,)*DF)*4.*DE SUMS = SUMS+S(I)*RADFAC SUMA = SUMA+A(I)*RADFAC LOOPP SUMB = SUMB+B(I.)*RADFAC NUMS = 2.*(S(N)-S(N-1))*N NUMA = 2o*(A(N)-A(N-1))*N NUMB = 2.*(B(N)-B(N-1))*N DENS = S(N)-SUMS DENA = A(N)-SUMA DENB = B(N)-SUMB NUSTED = NUMS/DENS PRINT COMMENT$OTHE STEADY NUSSELT NUMBER IS $ PRINT RESULTS NUSTED PRINT COMMENT$OTHE STEADY PERIODIC NUSSELT NUMBER IN PI/6 INC 1REMENTS IS $ SUMNUS = O0 THROUGH LOOPS,FOR ARG=PI/6C,PI/6.,ARGG.e2.*PI NUSPER = (NUMS+NUMA*SIN.(ARG)+NUMB*COS!(ARG))/(DENS+DENA l*SIN.(ARG)+DENB*COS.(ARG)) SUMNUS = SUMNUS+NUSPER LOOPS PRINT RESULTS-ARG.NUSPER TANUPR = SUMNUS/12. PRINT COMMENT$OTHE TIME AVERAGE OF THE PERIODIC NUSSELT NUMBE 1R IS $ PRINT RESULTS TANUPR PRINT RESULTS S(O-).O-S(-NI',V-(O)...V(M),PHIF(O).o.PHIF(N), 1PHIW(0)..oPHIW(M),AMPRAF(O)'..AMPRAF(N),AMPRAW(O)...AMPRAW(M) LOOPB CONTINUE TRANSFER TO START VECTOR VALUES TITLE = $1H1,S30',57HHEAT TRANSFER TO LAMINAR FL lOW IN TUBES WITH INTERNAL HEAT/S30,62HGENERATION - INCLUDING 2RADIAL DISTRIBUTION OF WALL TEMPERAtURE///*$ VECTOR VALUES XOUT=$1HO,20HCURRENT VALUE OF X = F10.6//*$ VECTOR VALUES ITTAB = $1HO,53HFOR THIS AXIAL LOCATION THE NUM 1BER OF ITERATIONS WAS I5//*$ END OF PROGRAM....

PROGRAM B - SMALL r DIMENSION A(25),B(25),E(25),F(25),G(25) H(25),S(25),U25), 1W(25),PH'F 25 )PERCOF (25), OPPF(25),ADJF (25), 2AMPRAF (25 ) INTEGER -T iJ,K,M,N,P,COUNT PRINT COMMENT $2$ PRINT FORMAT TITLE START READ DATA PRINT COMMENT $1NEW DATA$ PRINT RESULTS EPSQ,BETAGAM,LPDELTA DELT = 10.OE-36 ALPHA = *001 PI = 3.14159 THROUGH LOOPTFOR N = 20,-1,N.L. 10 WHENEVER GAM*(((DELTA.P.2.) +2.*DELTA)/BETA +1./N).G.4.*N-2 PRINT COMMENT$OTHE NUMBER OF NODAL POINTS IN THE FLUID IS$ PRINT RESULTS N TRANSFER TO ENTRYA OTHERWISE CONTINUE LOOPT END OF CONDITIONAL PRINT COMMENT$ODIVERGENCE OCCURS IN THF WALL EQUATION WHEN N lIS REDUCED TO 10. DATA IS REJECTED $ TRANSFER TO START ENTRYA THROUGH LOOPU, FOR P=P,5,P/L G.o 5.0E3 WHENEVER P/L.G. (GAM*(N.P.2.)/(2**N-1)) PRINT COMMENT$OTHE NUMBER OF AXIAL NODAL POINTS IS$ PRINT RESULTS P TRANSFER TO ENTRYB OTHERWISE CONTINUE LOOPU END OF CONDITIONAL PRINT COMMENT$ODIVERGENCE OCCURS IN THE F-LUIb- EQUATIONS WHEN 1P IS INCREASED TO P/L = 5000. DATA iS REJECTED.$ TRANSFER TO START ENTRYB CONTINUE CA = N.P.2. CB = P/L CC = CA/GAM CD = CB/G'.,M DE = 1./N DF = 1./(N.P.3.) DG = 2.*N-1. DH = (DELTA.P.2,)+2.*DELTA DI = DH/DG DJ = GAM*((DH/BETA)+DE) PRINT COMMENT$OTHE CONSTANTS GENERATED BY THIS DATA ARE$ PRINT RESULTS CA,CBCC,CD,DE,DF,DGDH,DI,DJ THROUGH LOOPA,FOR I=0,1,I.G.N A(I) = 0. B(I) = 0. S(I) = 0. E(I) = 0. F(I) = O. LOOPA U(I) = 0. THROUGH LOOPB,FOR J=1,1,J.G.P THROUGH LOOPC,FOR I=0,1,I.G.N G(I) = O. H(I) = O. LOOPC W(I) = 0.

-125COUNT = 0 — ENT1RT'-T-T —TDZ-JYA —TTT1 —--— G7DVT* A(- N-'') -D-F* EP SQ- AJ A(N) = -(DG/DJ)*B(N) +(DG/DJ)*B(N-1) SN) = S(N-1)+D! THROUGH LOOPFFOR I=N-1.-1,I.L.1 -----------— A —1- c —T?7T2;F-'A-+ {I +1) +CC* ( 1.-1. / ( 2 *I ) *A (I - 1 ) + lB(I)+CD*(1.-(I.Po2. )C)E)/CA*EI))/(CD*(l-(IP.2.)/CA)+2**CC) B(I) = (CC*(I.+1./(2.*I))*B(I+1 )+CC*(i.-1./(2.*I))*B( I-1)1A(I)+CD*('lo-(I.oP.2.)/A*F(I( ) )/(CD*(1.-(I.PA2s)/CA)+2o*CC). LOOPF —..STT - "-CCAiT-'; — -;' - CB* (. - ( I. P 2 ) / CA +2.*CA ) )S ( I+1 ) 1+(CA*(1.-lo/. 2 )/2*I))/(CB*(l.-(I.P2)/CA)+2.*CA))*S(I-1) -' —-`- +Y2-I-CB-;-'TP-2-;T)-C('i'-/'('-B-:cB';'-;(' I -.P ~ 2)/CA)+2~*CA ))*U(I) A(O) = (4.*CA*A(1)+GAM*B(O)+CB*,E(O))/(4o*CA+CB) - ~ ~~~~~BTCOl- 143 CN*BTFlT-GA]1IAYO) +~~BCF -0)) / 1 —4 *CA+CB) S(O) = (4.*CA/(4.*CA+CBH)*S(1)+(CB/(4.*CA+CB))*U(O).-..-R 0 — 1 --'-TROUG-WOURGthPUV;F 1:-7. —T- T;G N WHENEVER *ABS.(A(I)).L.DELT.OR.~ABS.(B(I)).LDELT.OR..ABS.(.............'IISt-) —tr.L.1;ELrT —-— T — -—...-. TRANSFER TO ENTRYE - -- ---- UTEIW'T-SE ------------------ WHENEVER.ABS'.((A(I)-G(I))-/A(I))*L.ALPHA*AND.*ABS.((B(I-)-H(I.....'..), /B'(I) ).L -ALPHA.ANDABS ( (S (I )-W( I ) )/S( I ) ).L.ALPHA CONTINUE'-' -— *-THER-WI SE COUNT = COUNT + 1 ~~~CRE —--- —............. —-------—.........l —----—. —- - PRINT COMMENT $ *******************************************$ PRINT COMMENT$ CONVERGENCE IS TOO SLOW. DATA IS REJECTED$ PRINT RESULTS S(O)o~~S(N),W(0)...W(N)}A(O):.*A(N),G(O)e.*G(N) 1,B(O).'.B(-N),H(0).*H(N) PRINT COMMENT $ *******************************************$ -T- r 1TRANFI~E RT o —S TA'R-T — OTHERWISE CONTINUE END OF CONDITIONAL THROUGH LOOPL,FOR K2O,1,K.GN.. G(K) = A(K) - H(K-) -- B('K).. —-- ---------... —----------------- LOOPL W(K) = S(K) tRANSFER TO ENTRYC END OF CONDITIONAL'-FNTRYE... —-—.ED — DO-F —CONDTTTONALLOOPJ CONTINUE..- R.O.... —-...TCroUGH FEoOPF-UFORT — ~;T. G5jSE(I) = A(I) F(I) = B(I) LOOPM U(I) = S(I) X = J/CB PRINT COMMENT $0$ PRINT FORMAT XOUTX PRINT FORfiAT ITTAB,COUNT THROUGH LOOPN-iFOR I='O,.I-I;.G.N NORMF = 1./(.ABS*(B(I)). ----------------— OPPFITT- — NOR-MF*-B —T1.ADJF(I) = NORMF*A(I)............PHIF IT --— ATNI —-. OPPFCIT'A-D-JF'I r)-)....'...... LOOPN PERCOF(I) = (SQRT.((OPPF(I).P.2.)+(ADJF(I).P.2)))/NORMF...T............-.-..HROUGtH'-:OtOPFR —-.'.I.T.'. LOOPQ AMPRAF(I) = PERCOF(I)/(S(I)*EPSQ)

-126SUMS = S(O)/(2.*CA) SUMA = A(O)/(2.*CA) SUMB = B(O)/(2.*CA) - R OUGN-'HU P P0 -F F'R- -I= T:-F'.-G;-' —-.-.. RADFAC = (I*DE-(I.P.3&)*DF)*4.*DE -SU —--------------- S —-— SITTRRAD FA —— c - - SUMA S= UMA + A(I)*RADFAC -LCOPP.. —-- - -- SUMB- =-BSU —B [ — -):*'RA-DFA — RADFAC NUMS = 2**(S(N)-S(N-1))*N NUMA-*-= — 2' —ATN- -A ( N-1 )T*-Y-4N - NUMB = 2.*(B(N)-B(N-1))*N DENS = S(N)-SUMS DENA = A(N)-SUMA ----------------- D.E.N.B.;-BT.-.U.B ---------- NUSTED = NUMS/DENS P T- -- RTCOtMENT$OT HE —TEY-NUSS ELT' NUMBER IS $ PRINT RESULTS NUSTED' —-------------— PRNT —COMENT$THE —T-EA-D-PER-O -NU SSELT- NUMB ER I N P I /6 INC 1REMENTS IS $ -—'- U US —------ Uus ---- THROUGH LOOPS,FOR ARG=PI/6.,PI/6.,ARG.G. 2.*PI t ISPER- = — ('UMS+'NUMA*STNST4'"ARG) +NUMB*COS (ARG) ) /DENS+DENA* 1SIN.(ARG)+DENB*COS.(ARG)) -— NUS MNU +N —----— ER- - LOOPS PRINT RESULTS ARGNUSPER ----------------— T-TANUPR- = - SUMtUS / 1. 2PRINT COMMENT$OTHE TIME AVERAGE OF THE PERIODIC NUSSELT NUMBE PRINT RESULTS TANUPR ---------— PRIN T —RSUZTS-S O N-.HI —-CS1TFTnt-~.H.PF N I,.AMPRAF ( 0 ) ~ * ~ 1AMPRAF(N) -COOPRB- -----— CONT INU E - - - --- - - --- TRANSFER TO START VECTOR VALUES TITLE=$S30,i 57HHEAT TRANSFER TO LAMINAR FLOW IN 1TUBES WITH INTERNAL HEAT/530,54HGENERATION - WITH THE WALL TE....... 2MPERATURE LUMPED-RADIALLY ///*$ VECTOR VALUES XOUT = $1HO,20HCURRENT VALUE OF X =F106//*$ -------------— VECTOR'-VALUE-S —'TTAB-$=TIHO- -53HFOR —THIS AXI'AL LOCATION THE NUMBE 1R OF ITERATIONS WAS I5//*$ END-O-F -PROGRAM.-.. —.

APPENDIX D THIN WALL TEMPERATURE DROP In the experimental study described in Chapter IV a thin-wall tube was used for the test section. The computer solution for a thinwall tube assumes a "radially lumped" wall where it is assumed that the wall temperature does not vary from inside to outside at any axial location. The object of the analysis in this Appendix is to develop a method which uses the experimental data taken (outside wall temperature) to predict the inside wall temperature in order to check the validity of assuming a lumped wall. The system considered is shown in Figure 40 where the outer wall is insulated and the temperature there is known. Any arbitrary condition may exist at the inner wall of the tube. Energy is generated in the wall at the uniform volumetric rate q"= 9qo(I0 +6Sn ) while the measured outside wall temperature is given (in dimensionless terms) by Got C s t+ *Sen(re- *) The energy equation, also in dimensionless form, for the wall, neglecting axial conduction and assuming constant physical properties, is A a a' r r r A % ( J )= (l/a S 1 (D-2) -127

-128The problem can be divided into steady and periodic parts by assuming that cC- Qszge(Lr~o) lpy}(0Cy )6 such that r,,a _r _r__ ar 3ar (fQ 2) = 0 (D-2) flr (f, = = Gr and hap art r -A A G %,S (Ro, o) = O (D-3) ar R e' = Integrating Equation (D-2) twice and evaluating the arbitrary constants with the two boundary condition gives L/ISw 4- 4 2A (LnL R? /R (Dk4) Equation (D-3) is solved by assuming the existence of a fictitious

-129similar problem with a Cos PG generation rate / a _ 2Q I r a4-/ 4 - Cosre A 6 arL r ar A 3QaR R 1 3), = 0 (D-5) /R, ~c>) =Yr 3r Cos(rs- r) Multiplying Equation (D-3) by i and adding to Equation (D-5) defines a complex variable Z = Q + irpw whose differential equation is I OZ:. _ a~ p Ie A is 0' rr r o a (r.'F 0= 0 (D-6) i e (/ p; }= e e Finally, since the generation rate is periodic in its real and imaginary parts and since the asymptotic solution is desired, it is assumed that the time dependency of Z is also periodic and Z can be represented by Z (r), = r, e (D-7) where 8 is also complex. Substituting (D-7) into (D-6) yields iZ'9 /e a'8 /9C Jt ah~ la — or rR' -L A Ts~~~:*'~~

-130Eliminating e gives r~ r er h A AA,@ 3 @ _ % RD- 0(D-8) The general solution for Equation (D-8) is ~- =Aor)4+ B-o/r) + Kr After applying the two boundary conditions it is found that L, %kr V,( Io iT + C.) (D-9) The denominator of (D-9) may be reduced considerably. From Reference 17 Thus eZ = +3etro _' eFr9(L ei~+ (j(LY~f)(V (Lya)i @(D-l0) (D-10) r /.1/X>> /s/;^t,'B; 6re

-131The desired solution for *pw is given by Writing each factor in Equation (D-10) in its complex form (a + ib), multiplying all factors and taking the imaginary part gives (/X = (( F: Pi) Cos(lre- )+ (F- S;n(T - (D-ll) __ r((p - ) OSTO - (F? +F )SLnr) - C COsrG where R = - KEI, CTBERI - KER, B'? - BEI, FKER?7 f — BER, KErEI7 PL=, KERR, BERT - KEZIBE - BEI, (KEL7 -BERIERT The functions BER y, BEI y, KER y, KEI y, BER1 y, BEI1 y, KER1 y and KEI1 y are respectively, the real and imaginary parts of the complex variables (See McLachlan(l7) for details of these functions.) Rearranging (D-ll) into Sin rG and Cos PG terms +ear = ( f) Cos. ~ > Cos n (D-2) V-: vTP) rCS m(/+z l~i(-2

-132or Yp =- P3 Sir)Z' 4 CosJ7 (D-13) This can be finally expressed as ~p~' \/ 8L~ F 4 Sin(re- ) (D-14) where = tn-l -P4 The wall temperature drop is given by QAk,,, _= -S ( 4- P&fl9 - >( s= n(re(i The steady part of ALw P from Equation (D-4), is 2 iC/sllr (IR~3 )n RL + (R.; )i (D-15) The periodic part of Atw, from Equation (D-12), is C/ipur = (x tSinFr Cos 4 - (:p_ Cos(e SinQ k ( Lc (P ( ) CosI- ( - p Sn (B P^)I- r CosrC

-133Combining Sin Ir and Cos r9 terms "] Ps, (C>os6re Or pi = -Pssnrcs- eosre (D-17) where D t an p' e It can be shown that Limpg O Re since Lim FP - - P. and ~Lem (..,

APPENDIX E THE EFFECT OF AXIAL CONDUCTION In Chapter II in the formulation of the problem it was assumed that axial conduction was negligible in comparison with the enthalpy flux in the fluid energy equation. Had axial conduction been included in the analysis, there would have been an added term on the right side of Equation (A-4) which is a Tf In terms of the dimensionless variables this contribution would be (R1 )2 magnitude of this term Re Pr 6x2 must be compared with the contribution of the enthalpy flux, which is (l-r2) xIn the problem formulation, when axial conduction was neglected, the dimensionless axial variable was chosen to be x = (x )(Re Pr)-lo When axial conduction effects are included, a more significant comparison can be made by specifying the axial variable to be (R)e In this case Ri the ratio of enthalpy flux to axial conduction would be /-re~P~)(,-~3 ~9 CY (I~e Pr)(/-r' ) &(Y/R) _ ______(E-l) Since the computer results are in terms of the right-hand side of Equation (E-l), this form of the ratio is used in the comparison. Goldstein(6) points out that for x > 0.25 the fluid temperature increases linearly with x so that is zero in this region. At smaller and smaller values of x the relative importance of axial conduction grows and becomes most important at the entrance. Also, at the entrance, the contribution of the enthalpy flux decreases as r - loO since the -134

-135 - fluid velocity goes to zero there. Thus, at the entrance and near the wall the ratio of enthalpy flux to axial conduction can be expected to fall. For example, with the present computer solution operating under the conditions B = 1.16, A 26.9, 8 0.126 and r = 0.95, the ratio of enthalpy flux to axial conduction is 2.56(ReoPr)2 x 10-5 and at r = 1.0 the ratio is zero. However, at r = 1.0 the ratio of radial conduction to axial conduction is 2.85(Re.Pr)2 x 10-5. This ratio evaluated at x = 0.002 is 3.2(Re.Pr)2 x 103 and continues to increase and approach infinity as the axial temperature gradient approaches a constant value for large x. In the wall, the same order of magnitude of the ratio of radial conduction to axial conduction at x = 0' and x = 0.002 exist at the fluid-wall interface but at r = Ro/Ri there is a very thin layer where - is small and axial conduction becomes quite predominant, even br2 for a large Peclet number, Re o Pr. From the foregoing evaluations it can be seen that for the effect of axial conduction to be insignificant in relation to the other terms in the energy equations the Peclet number must be greater than 2000. This is in contrast to the results of Schneider(24) who shows the minimum Peclet number to be l00 for slug flow. It must be pointed out, however, that inclusion of axial conduction in the fluid. and wall energy equations would change the resulting axial temperature gradients in such a way to increase the ratio in Equation (E-1l) and thus reduce the minimum Peclet number required for neglecting axial conductiono It can be concluded that the minimum Peclet number as given by Schneider (for slug flow) is too small for Hagen-Poiseuille flow in the entrance of the tube.

BIBLIOGRAPHY 1o Allen, Do No de G., Relaxation Methods. New York: McGraw-Hill, 1954. 2, Brown, G. M. "Heat or Mass Transfer in a Fluid in Laminar Flow in a Circular or Flat Conduit." J. AIChEo, 6, No. 2 (1960), 179. 30 Dennis, S. C. R. and Poots, G. "An Approximate Treatment of Forced Convection in Laminar Flow Between Parallel Plates " Appl. Sci. Res., A5 (1956), 453. 4. Dusinberre, G. M. "Calculations of Transient Temperatures in Pipes and Heat Exchangers by Numerical Methods." Trans. ASME, 76 (1954), 421. 5. Faddeeva, Vo No Computational Methods of Linear Algebra. New York: Dover, 1959. 6. Goldstein, S. Modern Developments in Fluid Dynamics. II, Oxford: Clarendon Press, 1938. 7. Graetz, L. Annalen d. Physik., 18 (1883), 79. 8. Graetz, L. Annalen d. Physik., 25 (1885), 337. 9. Jakob, M. Heat Transfer, I, New York: Wiley, 1949. 10, Kays, W. M. "Numerical Solutions for Laminar Flow Heat Transfer in Circular Tubes. " Trans. ASME, 77 (1955) 1256, 11. Koppel, L. B. and Smith, J. M. "Laminar Flow Heat Transfer for Variable Physical Properties." Journal of Heat Transfer, Trans. ASME, 84 (1962) 157o 12. Langhaar, H. L. "Steady Flow in the Transition Length of a Straight Tube." Journal of Applied Mechanics, 64 (1942), 55. 13. Leveque, M. A. Ann. Mines, 13 (1928), 201. 14. Lyche, B. C. and Bird, R. B. "The Graetz-Nusselt Problem for a Power-law non-Newtonian Fluid."' Chem. Eng. Sci., 6 (1956), 35. 15. Maslen, So H. On Fully Developed Channel Flow: Some Solutions and Limitations and Effects of Compressibility, Variable Properties, and Body Forces, NASA TR R-34, 1959. 16. MeAdams, W. H. Heat Transmission. New York: McGraw-Hill, 1954. -136

-1375 17. McLachlan, N. W. Bessel Functions for Engineers. Oxford: Clarendon Press, 1961. 18. Minneapolis-Honeywell Heiland Div., Catalog No. HC-301, Denver, Colo., 1958. 19. Organick, E. I. A Computer Primer for the MAD Language. CushingMallory, Inc., Ann Arbor, Michigan, 1962. 20. Perlmutter, M. and Siegel, R. "Unsteady Laminar Flow in a Duct with Unsteady Heat Addition." Journal of Heat Transfer, Trans. ASME, 83 (1961) 432. 21. Pigford, R. L. "Nonisothermal Flow and Heat Transfer Inside Vertical Tubes." Chem. Eng. Progress Symposium, 51, No. 17 (1955), 79. 22. Prins, J. A., Mulder, J. and Schenk, J. "Heat Transfer in Laminary Flow Between Parallel Plates." Appl. Sci. Res., A2 (1951), 431. 23. Schenk, J. and Van Laar, J. "Heat Transfer in Non-Newtonian Laminar Flow in Tubes." Appl. Sci. Res., A7 (1958), 449. 24. 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. Sellars, J. R., Tribus, M. and Klein, J. S. "Heat Transfer in a Round Tube or Flat Conduit - The Graetz Problem Extended." Trans. ASME, 78 (1956), 441. 26. Siegel, R. "Heat Transfer for Laminar Flow in Ducts with Arbitrary Time Variations in Wall Temperature." Journal of Applied Mechanics, Trans. ASME, 82 (1960), 241. 27. Siegel, R. "Laminar Convection in a Channel with Wall Heat Capacity and Wall Heating Variable with Axial Position and Time." To be published in Journal of Heat Transfer, ASME. 28. Siegel, R. and Perlmutter, M. "Heat Transfer for Pulsating Laminar Duct Flow." Journal of Heat Transfer, Trans. ASME, 84 (1962), 111. 29. Siegel, R., Sparrow, R. M. and Hallman, T. M. "Steady Laminar Heat Transfer in a Circular Tube with Prescribed Wall Heat Flux." Appl. Sci. Res., A7 (1957), 386. 30. Singh, S. N. "The Determination of Eigen-values of a Certain SturmLiouville Equation and Its Application to Heat Transfer." Appl. Sci. Res., A7, 1958.

310 Singh, So N. "Heat Transfer by Laminar Flow in a Cylindrical Tube." Applo Sci. Res., A7, 1958. 32, Tribus, Mo and Klein, J. S. "Forced Convection from Non-isothermal Surfaces. " Heat Transfer: A Symposium, University of Michigan Summer Conference, UMRI, Ann Arbor, Michigan, (1953), 211. 33. Van der Does de Bye, J. A. W. and Schenk, J. "Heat Transfer in Laminary Flow Between Parallel Plates." Appl. Sci. Res., 43 (1952), 308. 34, Whiteman, I. R. and Drake, W. B. "Heat Transfer to Flow in a Round Tube with Arbitrary Velocity Distribution." Trans. ASME, 80 (1958), 728. 35. Yang, W. J. The Dynamic Response of Heat Exchangers with Sinusoidal Time Dependent Internal Heat, Generation. Ph.D. Dissertation, The University of Michigan, Ann Arbor, Michigan, 1960. 36. Zeiberg, S. L. and Mueller, W. K. "Transient, Laminar, Combined Free and Forced Convection in a Duct." Journal of Heat Transfer, Trans. ASME, 84 (1962), 141. 37. Morse, P. M. and Feschbach, Ho Methods of Theoretical Physics, II, New York: McGraw-Hill (1953), 1092. 3 9015 02227 0949