ENGINEERING RESEARCH INSTITUTE DEPARTMENT OF AERONAUTICAL ENGINEERING UNIVERSITY OF MICHIGAN ANN ARBOR APPLICATION OF DIFFERENCE TECHNIQUES TO HEAT FLOW PROBLEMS USING THE ELECTRONIC DIFFERENTIAL ANALYZER by R. M. Howe Assistant Professor of Aeronautical Engineering University of Michigan Project 2115-3-T Office of Ordnance Research, U. S. Army Contract No. DA-20-018-ORD-21811 AIR-10 May, 1954

PREFACE This report continues the description of the application of electronic differential analyzers to partial differential equations of physics and engineering by using difference techniques. Earlier reports dealt with the lateral vibrations of beams, including both linear1' 13 14 and nonlinear damping terms. In this report we are concerned with both linear and nonlinear dynamic problems involving heat flow. In all cases time is preserved as a continuous variable; distance through the conducting medium is divided into stations and derivatives with respect to the spacial coordinates are approximated by finite differences. The resulting simultaneous ordinary differential equations are then solved by the electronic differential analyzer. Heat-flow problems treated include one and two-dimensional flow in rectangular media as well as flow in cylindrical and spherical media. In all cases theoretical accuracy of the difference method as a function of the number of stations is analyzed and example computer solutions are presented. Change of independent variable to improve accuracy and to solve flow in semi-infinite media is discussed with sample solutions. The heat equation for a medium having a conductivity proportional to its temperature is solved with the electronic differential analyzer and results are compared with a particular exact solution. It should be pointed out that the heat equation can be solved by the difference method with a passive network of resistors and capacitors. The active-circuit approach presented here has the advantage of versitility, ability to handle nonlinearities, and lowimpedance outputs. In addition, the extreme flexibility'in time scale of the electronic differential analyzer should allow a real-time simulation of heat-transfer systems when those systems are part of an automatic control process. Such a simulation could be extremely useful in designing and testing the control system. One should also note that all conclusions reached in this report concerning theoretical accuracy of the difference method in solving the heat equation n c readily be applied to the wave equation. This follows from the similarity between the two equations, the only difference being in the time derivative term; the heat equation i

involves a first-order partial time derivative while the wave equation involves a second-order partial time-derivative. Thus the normalmode shapes for both equations will be identical, while the normalmode frequencies for the wave equation will be the square root of the equivalent exponential decay constants for the heat equation. For this reason one can obtain from the curves presented in this report the errors in normal-mode shapes and frequencies as a function of the number of stations for the wave equation in Cartesian, cylindrical, and spherical coordinates. The computer solutions were obtained on the electronic differential analyzer installations in the Department of Aeronautical Engineering. The computers used included an 80-amplifier facility with components essentially similar to those discussed in another 3 report and a REAC * facility. Solutions were recorded on an xy plotting table. The reader of the previous reports may question the inclusion in this report of background material on the theory of electronic differential analyzers. This is done so that one not familiar with this type of computer can read the report successfully without recourse to previous reports. * Reeves Electronic Analog Computer, Reeves Instrument Corp., New York 28, New York ii

DISTRIBUTION LIST Office of Ordnance Research Commanding Officer Box CM, Duke Station Rock Island Arsenal Durham, North Carolina Rock Island, Ill. Office, Chief of Ordnance Commanding Officer Washington 25, D. C. Watertown Arsenal ATTN: ORDTB Watertown 72, Mass. Office, Chief of Ordnance Chief, Bureau of Ordnance (AD3) Washington 25, D. C. Department of the Navy ATTN: ORDTX-P Washington 25, D. C. Office, Chief of Ordnance Commander Washington 25, D. C. U. S. Naval Proving Ground ATTN: ORDTX (Technical Library) Dahlgren, Virginia Chief, Detroit Ordnance District Director 574 East Woodbridge Applied Physics Laboratory Detroit 31, Michigan Johns Hopkins University 8621 Georgia Avenue Commanding General Commanding Generl Silver Spring 19, Maryland Aberdeen Proving Ground, Md. ATTN: BRL Corona Laboratories National Bureau of Standards Commanding Officer Corona, California Frankford Arsenal Bridesburg Station U. S. Naval Ordnance Laboratory Philadelphia 37, Penna. White Oak, Silver Spring 19, Md. ATTN: Library Division Commanding Officer Picatinny Arsenal The Director Dover, New Jersey Naval Research Laboratory Washington 25, D. C. Commanding Officer ATTN: Code 2021 Redstone Arsenal Huntsville, Alabama iii

Commander, U. S. Naval Ordnance Technical Information Service Test Station, Inyokern P. 0. Box 62 China Lake, California Oak Ridge, Tennessee ATTN: Technical Library ATTN: Reference Branch Commanding General Commanding General Air University Research and Engineering Command Maxwell Air Force Base, Ala. Army Chemical Center, Maryland ATTN: Air University Library Technical Reports Library Commanding General SCEL, Evans Signal Corps Lab. Air Research and Development Command Belmar, New Jersey Baltimore, Maryland ATTN: RDR Office of the Chief Signal Officer Engineering and Technical Division Commanding General Engineering Control Branch Air Research and Development Command Room 2B273, The Pentagon Baltimore, Maryland Washington 25, D. C. ATTN: RDD ATTN: SIGGD Commanding General Commanding Officer Air Materiel Command Engineer Research and Development Wright- Patterson Air Force Base Laboratories Dayton 2, Ohio Fort Belvoir, Virginia ATTN: Flight Research Lab. F. N. Bubb, Chief Scientist Scientific Information Section Research Branch NAC for Aeronautics Research and Development Division 1724 F Street, Northwest Office, Assistant Chief of Staff, G-4 Washington 25, D. C. Department of the Army ATTN: Mr. E. B. Jackson, Chief Washington 25, D. C. Office of Aeronautical Intelligence Commanding Officer U. S. Atomic Energy Commission Signal Corps Engineering Laboratory Document Library Fort Monmouth, New Jersey 19th and Constitution Avenue ATTN: Director of Research Washington 25, D. C. ATTN: Director, Division of Research iv

Director National Bureau of Standards Washington 25, D. C. CADO U. B. Building Dayton 2, Ohio ATTN: DADO-SA Chief of Naval Research c/o Navy Research Section Library of Congress Washington 25, D. C. Jet Propulsion Laboratory California Institute of Technology 4800 Oak Grove Drive Pasadena 3, California Chief, Ordnance Development Division National Bureau of Standards Washington 25, D. C. v

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN TABLE OF CONTENTS Chapter Title Page PREFACE 1 INTRODUCTION 1 1. 1 Basic Equations for Heat Flow 1 1. 2 Equations for One-Dimensional Flow 2 1. 3 Solution by Separation of Variables 5 1.4 Replacement of Partial Derivatives by Finite Differences 8 1. 5 Principles of Operation of the Electronic Differential Analyzer 11 2 ONE-DIMENSIONAL HEAT FLOW 15 2. 1 Equations to be Solved 15 2. 2 Boundary Conditions 15 2. 3 Initial Conditions 16 2.4 Complete Set of Difference Equations for One-Dimensional Heat Flow 16 2. 5 Electronic Differential Analyzer Circuit for Solving One-Dimensional Heat Flow 17 2. 6 Theoretical Accuracy of the Difference Method 19 2. 7 Universal Curves for Obtaining Temperature Distributions for Arbitrary Initial Conditions 24 3 HEAT FLOW IN TWO AND THREE DIMENSIONS; CARTESIAN COORDINATES 34 3. 1 Heat Equation in Three Dimensions using Cartesian Coordinates 34 3. 2 Solution by Separation of Variables in a Homogeneous Rectangular Medium 35 3. 3 Theoretical Accuracy of the Difference Method 37 vi

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Chapter Page 3 3.4 Solution in a Two-Dimensional Homogeneous Medium 38 4 HEAT FLOW IN CYLINDERS 43 4. 1 The Heat Equation in Cylindrical Coordinates 43 4. 2 Solution by Separation of Variable 44 4. 3 Theoretical Accuracy of the Difference Method for Cylindrical Heat Flow 47 4. 4 Analyzer Solution of Axially-Symmetric Heat Flow with Initially Constant Temperature 50 5 HEAT FLOW IN SPHERES 62 5. 1 Heat Equation in Spherical Coordinates 62 5. 2 Difference Equation for SphericallySymmetric Heat Flow 63 5. 3 Solution by Separation of Variables in the Spherically-Symmetric Case 64 5. 4 Theoretical Accuracy of the Difference Method for Spherically-Symmetric Heat Flow 65 5. 5 Analyzer Solution of Spherically-Symmetric Heat Flow Starting with Unit Temperature 65 5. 6 An Alternative Method of Writing the Difference Equations for Spherically Symmetric Heat Flow 71 6 CHANGE OF VARIABLE TO IMPROVE ACCURACY 74 6. 1 Regrouping of Station Locations 74 6. 2 Solution of Temperature Distributions in a Uniform Slab 74 7 SOLUTION OF HEAT-FLOW PROBLEMS IN SEMI-INFINITE MEDIA 80 vii

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Chapter Page 7 7.1 Equation to be Solved 80 7. 2 Exact Theoretical Solution 82 7. 3 Solution by the Difference Method Using a Change of Spacial Variable 82 8 NONLINEAR HEAT-FLOW PROBLEMS 87 8. 1 Introduction 87 8. 2 Nonlinear Problem to be Solved 87 8. 3 Formulation of the Difference Equations and Analyzer Circuit 88 8. 4 Analyzer Solution for Uniform Initial Temperature Distribution Across the Slab 90 8. 5 Exact Particular Solution by Separation of Variables; Comparison with the Difference-Equation Solution 90 APPENDIX I - CALCULATION OF DECAY CONSTANTS AND MODE SHAPES FOR THE DIFFERENCE APPROXIMATION TO AXIALLY-SYMMETRIC HEAT FLOW 97 APPENDIX II - CALCULATION OF DECAY CONSTANTS AND MODE SHAPES FOR THE DIFFERENCE APPROXIMATION TO SPHERICALLY-SYMMETRIC HEAT FLOW 104 BIBLIOGRAPHY 109 viii

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN LIST OF FIGURES Figure Page 1-1 Temperature Distirubtion Between Two Slabs 3 1-2 Temperature Distributions Across Conducting Slab 9 1-3 Station Arrangement for N = 10- 1/2 10 1-4 Operational Amplifiers 12 1-5 Schematic of Servo Multiplier 14 2-1 Analyzer Circuit for One-Dimensional Heat Flow 18 2-2 Analyzer Circuit for Heat Flow Through a Homogeneous Medium 20 2-3 Analyzer Solution for Heat Flow in a Uniform Slab with Unit Initial Temperature Distribution 21 2-4 Comparison of Analyzer and Theoretical Solution to Heat-Flow Problem 22 2-5 Percentage Deviation in Normal-Mode Constant n for One-Dimensional Heat Flow 25 2-6 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at All Stations 26 2-7 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at Station 1 28 2-8 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at Station 2 29 2-9 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at Station 3 30 2-10 Temperature versus Time in a Uniform Slab with Unit InitialTemperature at Station 4 31 2-11 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at Station 5 32 2-12 Temperature versus Time in a Uniform Slab with Unit Initial Temperature at Station 6 33 3-1 Two-Dimensional Problem in Heat Flow 40 3-2 Analyzer Circuit for Two-Dimensional Heat Flow in an Isotropic Homogeneous Medium 41 ix

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN LIST OF FIGURES (Continued) Figure Page 3-3 Analyzer Solution for Two-Dimensional Heat Flow in a Square Isotropic Homogeneous Medium 42 4-1 Percentage Deviation of Decay Constant as a Function of Number of Cells for Radial-Dependent Modes of Cylindrical Heat Flow 51 4-2 Comparison of 4-1/2-Cell Mode Shapes with Jo, Modes 1 and 2 52 4-3 Comparison of 4-1/2-Cell Mode Shapes with Jo' Modes 3 and 4 53 4-4 Axially-Symmetric Temperature Distributions in a Uniform Cylinder Following Unit Initial Temperature Distribution 55 4-5 Axially-Symmetric Temperature Distribution with l/z (0) = 1 56 4-6 Axially-Symmetric Temperature Distribution with 1-1/2 (0) = 157 4-7 Axially-Symmetric Temperature Distribution with Uz-1/2 () 1 58 4-8 Axially-Symmetric Temperature Distribution with u3-1/2(0) = 1 59 4-9 Axially-Symmetric Temperature Distribution with U4-1/2 (0) 60 4-10 Axially-Symmetric Temperature Distribution with u5-1/2 (0) 1 61 5-1 Percentage Deviation of Decay-Constant for Spherically Symmetric Heat Flow as a Function of the Number of Cells 66 5-2 Comparison of 4-1/2-Cell and Continuous First and Second-Mode Shapes for Spherically Symmetric Heat Flow 67 x

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN LIST OF FIGURES (Continued) Figure Page 5-3 Comparison of 4-1/2-Cell and Continuous Third and Fourth-Mode Shapes for Spherically-Symmetric Heat Flow 68 5-4 Comparison of 10-1/2-Cell and Continuous First and Second-Mode Shapes for Spherically-Symmetric Heat Flow 69 5-5 Comparison of 10-1/2-Cell and Continuous Third and Fourth- Mode Shapes for Spherically-Symmetric Heat Flow 70 5-6 Analyzer Solution of Spherically-Symmetric Heat Flow with Unit Initial Temperature and Walls at Zero Temperature 72 6-1 Comparison of Station Locations for the Space-Variable Transformation y = 4 75 6-2 Analyzer Solution for Heat Flow in a Uniform Slab; Equal Stations along y, where y = x 77 6-3 Comparison of Analyzer and Theoretical Temperature Distributions Across Uniform Slab with Unit Initial Temperature 79 7-1 Semi-Infinite Bar 80 7-2 Temperature Distributions along a Semi-Infinite Bar at Various Times 83 7-3 Analyzer Solution for Temperatures in a Semi-Infinite Bar 85 8-1 Analyzer Circuit for the Temperature at the nth Station; Conductivity Proportional to the Temperature 89 8-2 Analyzer Solution for Heat Flow with Conductivity Proportional to Temperature; Conducting Wall Held at Temperature U = 0. 91 8-3 Analyzer Solution for Heat Flow with Conductivity Proportional to Temperature; Conducting Wall Held at Temperature Uo = 0. 2 92 xi

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN LIST OF FIGURES (Continued) Figure Page 8-4 Analyzer Solution for Heat Flow with Conductivity Proportional to Temperature; Conducting Wall Held at Temperature U = 0. 5 93 8-5 Analyzer Solution X (x) for the Nonlinear Eigenvalue Problem 95 8-6 Analyzer Solution for Checking the Exact Solution by Separation of Variables 96 xii

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 1 INTRODUCTION It seems hardly necessary to point out the importance of the heat equation and the many problems in heat transfer to which it is applicable. Since the heat equation is a partial differential equation containing derivatives with respect to both time and spacial coordinates, it is necessary to convert it to one or more ordinary differential equations before solution by means of the electronic differential analyzer. This is because the electronic differential analyzer can integrate with respect to only one variable, namely time. If the original partial differential equation is linear and has suitable boundary conditions, it can be converted to ordinary differential equations by separation of variables. The resulting eigenvalue equations can be solved directly by the 4-9 electronic differential analyzer. If the heat-flow equation is nonlinear, the technique of separation of variables cannot be used due to the failure of the superposition principle. However, by considering the temperature only at discrete stations within the conducting medium and by approximating spacial derivatives with finite differences, the original nonlinear partial differential equations can be converted to a system of nonlinear ordinary differential equations capable of being handled by the electronic differential analyzer. Even when the problem is linear, the difference technique, although requiring more equipment, gives quicker and more direct answers and easily allows the introduction of time-varying boundary conditions, nonuniform thermal characteristics through the medium, arbitrary time-dependent heat sources, etc. The directly-recorded voltage outputs of the electronic differential analyzer represent the temperature and heat flux at each station as a function of time. 1.1 Basic Equations for Heat Flow The basic equation of heat flow is given by c~au (1-1) cu = V KVu +S at 1 -

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where u = temperature and is a function of the spacial coordinates and time, K = thermal conductivity, in general, a function of the spacial coordinates, c = specific heat, a function of spacial coordinates, 6 = density, also a function of spacial coordinates, t = time, S = rate of heat supplied per unit volume by sources, in the medium, a function of spacial coordinates and time. The left-hand side of equation (1-1) represents the rate at which heat is stored in a unit volume due to the heat capacity of the medium. The right-hand side represents the rate at which the unit volume receives heat, first due to heat conduction into the volume from the neighboring medium (the V * K V u term) and second, due to the heat flow into the volume from sources within the volume itself (the S term). The conductivity times the gradient of the temperature (- K Vu) is a vector representing the heat flux. The components of - K Vu represent the heat flow through a unit surface normal to the direction along which the component is taken. In a given heat-flow problem it is necessary to stipulate spacial boundary conditions either on the temperature u or the heat flux - K V u, as well as the initial temperature distribution throughout the medium. 1. 2 Equations for One-Dimensional Flow In order to simplify the discussion to follow, let us assume that spacial variations in the temperature u are confined to a single direction along which the coordinate is x. Equation (1-1) then becomes (x) au(x,t) = a [K(x) U(xt)] +s(x, t) (1-2) at ax ax where we have defined the heat capacity C (x) by C (x) = c () 6 (x) (1-3) For example, u (x, t) in Equation (1-2) could represent the temperature distribution in a medium between two infinite slabs, as shown in Figure 1-1. Let us __________________________-_ 2

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN TEMPERATURE u(x t) U UQ \NSULATOR =,CONSTANT X=O x=L Figure 1-1. Temperature Distribution Between Two Slabs. assume that the medium is bounded at x = 0 and at x = L. The heat flux F- along the x direction is given by au F- = - K (1-4) x It is convenient to define a dimensionless distance variable x such that the distance through the conducting medium is unity. Thus let x = L (1-5) from which a a dx 1 a a_. = _dx = _ (1-6) ax ax dx L ax Also let us transform the variable characteristics of heat capacity C (x) and conductivity K (x) into dimensionless variables 0C (x) and K (x) respectively. Thus let ---------------------------- 3 --------------------------— ~~~ ~~ ~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN C (Lx) = CO C (x) (1-7) and K(Lx) = K0 K (x) (1-8) Here Co and K are constants equal to the maximum value of C (x) and K (x) respectively. For a homogeneous medium 9C (x) = fK (x) = 1. From Equations (1-6), (1-7), and (1-8) the one-dimensional heat equation (1-2) becomes CO au(x, t) 1 a au(x,t) 1 0C (x) - = - K (x) + —S (Lx,t) (1-9) K at L ax ax K O O Next we introduce a dimensionless time variable t given by K t = t (1-10) C L o from which a a dt K 8 a - - - - - -Z. - (1-11) at at dt C0L at In terms of t Equation (1-9) becomes C ( au(xt) a au(x,t) + S (x,t) ~c (x) -a — x K (X) - - +S(Xt ( ) where 2 C L L O S (x,t) = S(Lx, K t) (1-13) O O For the special case of a uniform medium with no internal heat sources we have the familiar equation au = a2u (1-14) at ax Equation (1-12) must have boundary conditions specified. For example, if at all times the left side of the medium in Figure 1-1 is held at zero temperature while the right side is insulated (no heat flow past the wall), the boundary conditions become u (o,t) = 0 (1-15) and au(l,t) = 1-16) K- (1) -ax

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Let us denote the initial temperature distribution at t = 0 by u (x,o) = U(x) (1-17) 1. 3 Solution by Separation of Variables Before considering the solution of Equation (1-12) by the difference method let us review the solution by separation of variables. This will not only give us complete theoretical solutions for specific problems against which we can later check solutions by difference methods, but it will also allow us to compare the normal modes for the continuous solution with those using the difference technique. This latter comparison will, in turn, allow us to estimate the accuracy of the difference method for more general problems. Consider the case where S (x,t) = 0. Assume that the temperature u (x,t) can be written as u (x,t) = X (x) T (t) (1-18) Substituting Equation (1-18) into Equation (1 - 12) we have (x) (x) dT(t) d dX(x) (1-19) dt d-x-K(x) dx or 1 dT(t) 1 d dX(x) (1-20) T(t) dt 0c(x)X(x) dx dx Since the left side of Equation (1-20) is a function only of t and the right side is a function only of x, the only way they can be equal for all x and t is for both to be equal to the same constant. Thus let 1 dT(t) _ T(t) dt or dTt +3T (t) = 0 (1-21) In the same way d d(x) + p c(x) X(x) = 0 (1-22) dRx K(x) dx C By separating the variables we have transformed the original partial differential equation (1-12) into two ordinary differential equations (1-21) and (1-22). The solution of Equation (1-21) is simply an exponential function with decay-constant (, i. e., T (t) = Ae-Pt (1-23) ---------------------- 5 ---------

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where A is an arbitrary constant. Equation (1-22) is subject to the boundary conditions of the original problem, since homogeneous boundary conditions on u (x,t) must also be met by X (x). For example, for the medium with zero temperature at x = 0 and zero heat flux at x = 1 the conditions on X (x) are from Equations (1-15) and (1-16) X(0) = K dX(1) = 0 (1-24) 0) O ~K("dx The solutions to Equation (1-22) can satisfy the boundary conditions (1-24) only for discrete values n of the constant. The solution X1 (x) corresponding to the smallest allowable value p, gives the temperature distribution through the medium corresponding to the lowest mode of exponential decay with time. X2 (x) gives the shape of the second mode, which decays exponentially with decay-constant P2, etc. The discrete values of Pn are known as eigenvalues, and the corresponding solutions Xn (x) are known as eigenfunctions or normal modes. The higher the mode (i.e., the larger the integer n), the larger the decay constant n, and the faster will that particular mode shape X (x) decay or die out. One can show that any arbitrary temperature u (x, t) can be represented by the proper combination of the normal modes. Thus 00 u (x,t) =- A X (x)e en (1-25) (' nil The constants A can be evaluated by applying the initial condition of Equati (1-17) and by use of the orthogonality properties of X (x). Equation (1-22) is known as the "Sturm-Lionville Equation, " and it is easy to show from the equation itself and the boundary conditions that10 f (xC x)X(x) (x) (x) dx = 0 n m (1-26) N n = m n where n and m are integers. At t = 0 Equation (1-25) becomes 00 U (x) = AnX (x) (1-27) n=1 n n Multiplying both sides of Equation (1-27) by C- (x) Xm (x) and integrating with respect to x from 0 to 1 we have __________________________________ 6 _________________________________

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 1 oo 1 Sf 0C (x ) U () (x) dx - An f!C (x) Xn (x) X (x) dx (1-28) 0f ~c~x)U~x)X n=1 0 From the orthogonality relationship in Equation (1-26) it is clear that all the integrals on the right side of Equation (1-25) vanish except the one where m = n, from which the right side is A Nn. The solution for A is then A = (x) U (x) X (x) dx (1-29) n Nn 0 n nO A method for calculating the complete solution to the heat equation, at least when no internal heat sources are present, is now clear. First the normal-mode decay constants An and the corresponding normalmode functions X (x) must be computed from Equation (1-22). Next the mode coefficients A are found from Equation (1-29). Finally, the complete solution is obtained by summing the contributions of all the normal modes, as indicated in Equation (1-25). Although in theory an infinite number of these modes is required to represent exactly the initial temperature distribution U (x), in practice only the first few modes are usually required to give a suitably accurate representation for a reasonably well-behaved U (x) function. The electronic differential analyzer can be used to solve the eigenvalue equation (1-22) for the mode shapes and eigenvalues. To do this we let time on the analyzer represent distance x through the medium. In the same way the analyzer can be used to calculate the coefficients An from Equation (1- 29). Let us now consider the simple case where the medium is homogeneous and therefore 0K = 0C = 1. Equation (1-22) becomes 2 dX dx2 + pX = 0 (1-30) dx which has sin- x or cos p x as a solution. For the boundary conditions of Equation (1-24) to be met it is clear that the solution is X (x) = sin (n - 1/2) Tx, n = 1, 2, 3,.... (1-31) from which pn = (n - 1/Z)2 2 (1-32)

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN From the Xn (x) functions and pn values given above the complete solution can be written as the infinite series given in Equation (1-25). As a specific problem consider the case where the initial temperature throughout the medium is equal to a constant Uo. Then from Equations (1-26), (1-29), and (1-30) we can solve for the coefficients An, obtaining 2U A = (1-33) n ir(2n-l) Thus the complete solution is given by 2U2 u(x, t) = o E sin(n-l/2)wrx -(n-l/2)t (1-34) uT n=1 (2n-1) This solution actually represent an infinite number of sinusoidal temperature distributions across the medium from x = 0 to x = 1. At t = 0 the sine waves all add up to give the initial flat distribution. For t > 0 the sine waves decay exponentially at different rates, with the decay rates faster for those sine waves having more nodes and loops. The resulting temperature distributions at various times are shown in Figure 1-2. Note that after a short time has elapsed the distribution is almost pure sinusoidal representing the first mode. This is because the higher modes decay so rapidly. Later we will compare electronic differential analyzer solutions of these same problems using difference techniques with the results shown in Figure 1-2. 1. 4 Replacement of Partial Derivatives by Finite Differences In the one-dimensional heat flow problem considered in the previous section the temperature a was a function both of time t and distance x across the medium. Instead of measuring the temperature u at all distances x, let us measure u only at certain stations along x, as shown in Figure 1-3. Let ul be the value of u at the first x station, u2 be the value of u at the second x station, un be the value of u at the nth station. Further, let the distance between stations be a constant Ax. Thus u (x, t) is replaced by ul (t), u2 (t),.... etc., and we can approximate OK au/ax at the n-1/2 station by au n-1/2 auK I (u- n-lu ) (1-35) n-1/2 x 8 __________________________

1.0 0.002 0.04 0.80.16 0.6 0.53 2 U U0 0.4 0.0 0.0.4t 0.6 8 1. 0.2 0.0 I I I I 0.0 0,2 0.4 0,6 0,8 1.0 DISTANCE X THROUGH THE SLAB Figure 1-2. Temperature Distributions Across Conducting Slab. 9

- ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN U3U4 Us /? / \\Ul' INSULATOR U U0 / - O: CONSTANT Figure 1-3. Station Arrangement for N = 10-1/2. Indeed the limit of Equation (1-35) as Ax-.0 is just the definition of the partial x derivative at that point. In the same way we can approximate (a/ax) (K au/ax) at the nth station as a i U-KaJ au L-nl/ Kau n (1-36) ax OK ax n - A-X LK a n+l/2 K ax In-/ Thus the equation of heat flow balance at the nth station becomes from (1-12), (1-35), and (1-36) du 1 rKn/ C - n = F ['Wz (Un+l -n) - /K (un -u n) + f (t) 10.... ----------

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN If the distance between x = 0 and x = L is divided into N segments of length Ax (for example, N = 10-1/2 in Figure 1-3), then Ax = (1-38) Introducing a new time variable T given by T = Nt = 1 t, (1-39) (Ax)Z we have from Equation (1-37) du C dn dT = i/ (Un+l Un) - n-l/2(u -u ) +u "(T) (1-40) where n (T) = S T) () -41) Equation (1-40) is iterated for different values of n until the boundaries at x = 0 and x = 1 are reached. Boundary conditions will be discussed in the next chapter. Note that the set of simultaneous first-order equations represented by (1-40) at each station consist of ordinary differential equations; the only derivatives are with respect to the time variable T. Thus the equations can be handled by the electronic differential analyzer. 1. 5 Principles of Operation of the Electronic Differential Analyzer Although in previous reports in this series we have reviewed briefly the principles of operation of the electronic differential analyzer, to lend continuity to the present report this material is again included in this section. The reader who wishes to become more familiar with this 411 type of computer is directed to other references.' Those already familiar with the electronic differential analyzer may omit this section. The basic unit of the electronic differential analyzer is the operational amplifier, which consists of a high-gain do amplifier having a feedback impedance Zf and one or more input impedances, as shown in Figure 1-4. To a high degree of approximation the output voltage eo of an operational amplifier is equal to the input voltage multiplied by the ratio of feedback to input impedance, with a reversal of sign (Figure 1-4a). If several input resistors are used, the output voltage is proportional to the sum of the input voltages (Figure 1-4b). If an input resistor and feedback 11

-i a.) OPERATIONAL AMPLIFIER Ra ia Rf if ea Rb lbl OUTPUT eb I e2 Rc ic ec 2 = (Raa Rbb R c) b.) OPERATIONAL AMPLIFIER AS A SUMMER i2 C INPUT R ilt; o OUTPUT 1 e2 e2 -- RC e, dt C.) OPERATIONAL AMPLIFIER AS AN INTEGRATOR Figure 1-4. Operational Amplifiers. 12

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN capacitor are used, the output voltage is proportional to the time integral of the input voltage (Figure 1-4c). The operational amplifiers shown in Figure 1-4 can therefore be used to multiply a voltage by a constant factor, invert signs, sum voltages, and integrate a voltage with respect to time. To multiply several voltages a servomechanism which drives potentiometers is the most commonly used device. In Figure 1-5 the block diagram of a servo multiplier is shown. It consists of a number of linear potentiometers ganged together and driven by a servo motor. The reference voltage + VR is connected across one of the pots, and the variable tap voltage VR is subtracted from the voltage Z. The resulting error signal E = Z - aVR is sent through a high-gain servo amplifier and applied to the servo motor. The motor drives the variable tap in the proper direction to reduce the error to zero, i. e., to make aVR = Z. In this way the tap position on all of the ganged pots is proportional to the voltage Z. If ~ X and + Y are applied across each of the remaining two pots shown in Figure 1-4, it is apparent that the variable tap voltages will be XZ/VR and YZ/VR respectively. Thus the servo multiplier can generate output voltages proportional to the product of input voltages. For the electronic differential analyzer solutions obtained in this report REAC* Servo Unit S-101 Mod 4 servos were used. Accuracy of multiplication is about 0.1% of full scale (+ 100 volts). By employing operational amplifiers for summation and integration, and servos for multiplication, we are able to solve nonlinear heat-flow problems. * Reeves Electronic Analog Computer, Reeves Instrument Corp., New York 28, New York 13

X o <^z~~~-. F —-- VR Yc —-1 l -Xo,~ oI -YO I vR SERVO HIGH AIN R MOTOR SERVO --- AMPLIFIER VR o REFERENCE < POT -VR Figure 1-5. Schematic of Servo Multiplier. 14

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 2 ONE-DIMENSIONAL HEAT FLOW 2. 1 Equations to be Solved In this chapter we will solve by difference methods the onedimensional heat flow problem discussed in Sections 1. 3 and 1. 4 of the previous chapter. The equation to be solved is OC (x) aux, t) = (x) au(xt) + s (xt) (1-12) ac (x) uxt0 - ~ 0K (x) ax with boundary conditions at x = 0 and x = 1 on the temperature u or the heat flux OK au/ax. As discussed in Section 1.4, we can convert the above partial differential equation into a system of simultaneous ordinary differential equations by considering the temperature u only at equally spaced stations along x. The equation describing heat balance at the nth station is du w C dT = K (Un+l - n) - OK /(un - n-l + A) n (T) (1-40) n n+l/2 n-l/Z where T = N2t, N being the number of stations or cells and where n () = (1/Nz) S ( ). N 2. 2 Boundary Conditions Let us now consider the representation of the boundary conditions when using the difference method. If at x = 0 the temperature uo is specified, then this boundary condition is imposed by letting u0 equal the specified value in the difference equation for ul [n = 1 in Equation (1-40)] Note that the boundary condition at x = 0 may require that uo equal zero, a fixed constant, or a known time-varying function. Thus the problem of timedependent boundary conditions is readily handled. If at x = 1 the boundary condition specifies the temperature, then uN, the prescribed temperature at x = 1, is introduced into the difference equation for N_1 [n = N-l in Equation (1-40)]. On the other hand, the heat flux may be specified at x = 0 or x = 1 as a boundary condition. We recall from Equations (1-4), (1-5), and (1-8) that 15

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN au K 8u F- = -K (x) - - K (x). (2-1) Xx L ax where F- is the heat flux in the x (and hence x) direction. Approximating x au/ax by a finite difference we have K N F- / K (Un+1- U) (2-2) X n+1 /2 L n+1/2 as the expression for the heat flux at the n+l/2 station. If the flux must be equal to a specified value, say FN at x = 1, then from Equation (2-2) L K (UN+l/2 - uN-/2) f = FN N K N This expression for OKN (uN+1/2 - uN-1/2) = fN is used in the difference equation for UN 1/2 [n N-l/2 in Equation (1-40)]. Clearly the required flux FN at the boundary may be zero, a fixed constant, or a known timevarying function. Note that if we consider the temperatures as occurring at integral stations, the heat flux occurs at half-integral stations. If the boundary condition is on the heat flux, the boundary therefore occurs at a half-integral station. On the other hand, if the temperature is specified at the boundary, the boundary occurs at an integral station. Thus in Figure 1-3 the left-hand boundary, at which uo is a constant, occurs at the zero station, while the right-hand boundary, at which the flux is zero, occurs at the 10 + 1/2 station. Hence N = 10-1/2 in Figure 1-3. 2. 3 Initial Conditions The initial temperature distribution at zero time was specified in Equation (1-17) as u (x,0) = U(x) (1-17) In the difference-equation formulation this becomes an initial temperature condition at each of the stations. Thus in general let u (0) = U = U(n/N) (2-4) n n 2. 4 Complete Set of Difference Equations for One-Dimensional Heat Flow Let us now write down the complete set of differential equations representing one-dimensional heat flow when the temperature is specified as uo at x = 0 and the heat flux is specified as fN at x = 1. From Equation 16

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN (1-40) and the boundary conditions discussed in Section 2. 2 we have du1 C -1 = 1 u2 ul] K K [ut - uO(T)] +0 (T) 1 dT -1/2 l/ (U - K (u - u) - (u - u)+ (T) C2 dT K-1/2 3 2) KI- 2 1/) 2 (2-5) duN-1-1/2 ~C l — l/d = KNi (uN-1/2 - UN —1/2) CN-1-1/2 dT N-1 - KN (UN-1-1/2 - N-2-1/2) + ON-1-1/2 (T) du dUN-1/2 CN-1 d = fN - KN (uN-1/2 - N-l-1/2) +0N-1/2(T) N-1/2 dT N-i The initial conditions require that ul (0) = U1, u2 (0) = U etc. 2. 5 Electronic Differential Analyzer Circuit for Solving One-Dimensional Heat Flow The electronic differential analyzer circuit for solving the simultaneous first-order equations given in (2-5) is shown in Figure 2-1. Note that the output of each successive row of amplifiers is reversed. This allows the necessary differences to be taken without the use of sign-reversing amplifiers. Note also that the quantity fn, which is related by Equation (2-3) to the heat flux Fn, is available at any half-station. Single resistors in the circuit represent the heat-capacity parameter OCN and the conductivity parameter OK at each station. The heat source variable 5 n at each station is introduced as a time-varying input voltage. Note that the integrator time scale is RC seconds in the circuit in Figure 2-1. This means that RC seconds of real time on the computer equals one unit of T. It is possible to reduce the required number of amplifiers per station from three to one. In many ways, however, the circuit of Figure 17

-_, (T) 1/2 _ *,(r) c \ /'(f3/2 - fl/2)-, K 3/2 R~CC'>t f V (f5/2 - f3/2) + (2 K5/2 yl)f3 (-r) u3 \ / (f5/2 3 DKN-2 I-1''^a'^WN-r. N-1/2 / ( fN f-I ) ^N-3/2 NI c I Figure 2-1. Analyzer Circuit for One-Dimensional Heat Flow. 18 I~~~~~DK

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 2-1 is simpler despite the increased number of amplifiers, since it does have the one-to-one correspondence between resistor values in the circuit and parameters in the problem. For heat flow through a homogeneous medium, however, a circuit using but one amplifier per cell is quite useful. In this case Cn = 0Kn = 1, and the circuit shown in Figure 2-2 is utilized. As an example solution, consider the heat-flow problem discussed earlier in Section 1. 2. Here a uniform slab of unit thickness has initially a unit temperature distribution throughout. The edge at x = 0 is held at zero temperature while at x = 1 the heat flux is zero, corresponding to an insulated boundary. The theoretical temperature distributions across the slab at various times are shown in Figure 1-2. This problem was set up by the difference method for N = 6-1/2 on the electronic differential analyzer. A recording of the temperatures at each of the 6 stations as a function of time is shown in Figure 2-3. For this solution the circuit of Figure 2-1 was used, with 9Cn = OKn = 1 and n = 0. An integrator time scale of 5 seconds was used for the solution, so that 5 seconds of time in the computer solution was equivalent to unit T, which in turn equals 1/(6. 5)2 units of the dimensionless time variable t. In Figure 2-4, points from the electronic differential analyzer solution of Figure 2-3 are plotted on the theoretical temperature distributions shown previously in Figure 1-2. 2. 6 Theoretical Accuracy of the Difference Method Now that we have seen how the partial differential equation of heat flow can be converted to a set of first-order ordinary differential equations by difference methods and how the electronic differential analyzer can be used to solve these equations directly, let us examine the accuracy of the approximate solutions obtained by the difference equations. We recall that the gradient of the temperature half-way between two stations was approximated by (l/Ax) times the difference between the temperatures at the two stations. One can show by means of a Taylor series expansionl2 that the error in the gradient introduced by this approximation varies as 1/(Ax)2 where Ax is the interval between stations. Thus we should expect the accuracy of the difference method to improve as the square of the numbe of stations. 19

-Uo(T) 0.5R \ 0.5R C r-AA^^f> -U2 R\ I Ul R O.5R C - U2 D-32(T) -' / 0.5R C X3(T) 0.5R R C R -UN-3/2 R Figure 2-2. Analyzer Circuit for Heat Flow Through a Homogeneous Medium. N-3/2O20 R R N —/2 20

C...::.:::j i'i:::: i-:::...^....: -............., —-..................,..~~~~~~~~'.., (... I 6 ~ ~ ~ ~ ~ ~ ~~~~~~~f J' CJ E 00 (1I". ENE DIETZGEN COO. 346 D 0 0.0 0 0 0 0 0 o 33 3IN.LVWt 1 _ lI~3/1).... Sp f++ g 0. l g; 0 T 4 q0 W.. g f.. T lit rsf..1

1.0 0.0025 0.01 0,8 0.16 0,6 U, U0 0.4 t=.64 ANALIZER SOLUTION 0.0 0.0 0.2 0.4 0.6 0.8 1.0 DISTANCE X THROUGH THE SLAB Figure 2-4. Comparison of Analyzer and Theoretical Solution to Heat-Flow Problem. 22

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Since one would like to use as few operational amplifiers as necessary, it becomes important to estimate the exact accuracy of the difference method when relatively few stations are used. We will do this by finding the normal mode shapes and decay constants for the differenceequation representation of the uniform slab. These will then be compared with the sinusoidal mode shapes and the decay constants given in Equations (1-31) and (1-32). For a uniform conducting slab the difference equation at the ith station becomes from Equation (1-40) du. = u. - 2uf +u 6) dT 1i+l i-l (2-6) dT where we have assumed no heat sources within the slab ( = 0). To find the normal-mode solutions we assume that the temperature ui varies with time as a. e T, where a. and X are constants. If this is true, then (2-6) 1 1 becomes a set of N-1/2 simultaneous algebraic equations for the case of zero temperature at x = 0 and zero heat flux at x = 1. The only nontrivial solution of this set of equations is obtained when the determinant of the coefficients vanishes. This determinant, when expanded, becomes a polynomial in X of order N-1/2. The polynomial will have N-1/2 positive real roots n, which are the decay constants for the N-cell systems. To solve this determinant for a specific N is tedious, and to solve it in general would be next to impossible. Fortunately the roots Xn can be found directly by the following procedure. Assume that the spacial mode shapes for the difference equations are the same as for the continuous equation, i. e., sinusoidal. If this is true, then for the temperature u. at the ith station we have by analogy with Equation (1-31) (n-./2)r -\ T u. = a sin (n-l/)ri e n = 1,2,3,.... (2-7) Substituting Equation (2-7) into Equation (2-b) and noting that u+u +. = 2a sin (n-l/2)-i os (n-1/2)w (2-8) -i+l 1+ ui-1 N N we have -X a sin (n-1/Z)rri -X T (n-1/2i enT(c n2)w e a sin e (2cos -2) 23

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN from which X = 2[ 1-cos (n-/2)i) (2-9) n N Since the time variable T = N' t, the decay constant Pn referred to the dimensionless time variable t is N X Thus n ^n = 2N2 (l-cos (n4/) (2-10) After expansion of the cosine function in a power series we have 22 2 n 22 n 12N In the limit of infinitely many cells the 3n values given in the above equation reduce to those given earlier in Equation (1-32) for the continuous solution. As we had predicted, the error in the decay constant In for the difference method varies as 1/N, where N is the number of stations. In Figure 2-5 the percentage deviation in decay-constant due to the difference method as a function of the number of cells N is shown. To summarize, we see that when the spacial derivatives of the one-dimensional heat equation in a uniform medium are replaced by finite differences, the resulting mode shapes agree exactly with those for the continuous solution, whereas the decay constants (eigenvalues) for each mode are somewhat smaller. This means that the higher modes will decay somewhat more slowly when the difference approximation is used. Fortunately the higher modes are usually of less importance than the lowest modes in most problems. 2. 7 Universal Curves for Obtaining Temperature Distributions for Arbitrary Initial Conditions We have already seen in Figure 2-4 that the electronic differential analyzer solutions of temperature distributions within a uniform slab are fairly accurate, even when only 6-1/2 stations are used. The case treated was for zero temperature at one side of the slab, zero heat flux at the other side, and an initial unit temperature distribution across the slab. In Figure 2-6 the solution of Figure 2-3 has been recorded with an expanded time scale. In Figure 2-7 the temperatures at each of the stations are recorded as a function of time for intial conditions u1 (0) = 1, u2 (0) = u3 (0) = U4 (0) = U5 (0) = u6 (0) = 0. Figure 2-8 gives the solutions 24

'*MOTA ^EaH euoTsuauITaI-auo JOJ Ud JuI3suo3 apolJ -TBetUloN uT UOITBIAaJ aS'Euaoad's-Z ainxrTI S'133 JO 0 38 3WlnN 001 08 09 OS 0o 0~ O0 01 8 9 S ~ 2 I 0'001 - 0'08 - 0'09 - O'Og - 0-0 - - o' oI/'/ / /O' l - 0// - 0 /c o 9,:, /// / / o 0 --- /- / / - / _ __ —. -,/Q) ~ ~ ~ ~ ~ ~ 9 /A/v,7A 0t/0'1 - //0 -

[.0 (.O. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~. 1~~~......... *:::: ^ ^ ^ ^ ^ ^ ^ ^ ^.^ ^ ^ ^^ ^i^ ^ ^ ^ ^ ^ ^ ^ ^ ^ i^ ^ -S m ^-y- ^ ^ ^ ^~~~~~~~~~~~~~~~~~~~~~~~~~:: 0%6~.6 -^r —"-^ ^ ^^1..1"^"^ — ^^^!!^ I[~~~~~~~~~~~~. 8'"' 0. 4 ~~~~~~~~~~~~~~~~~~~~~..............:0::Ai! i!! t~ ~~~~~~~~~~~~~~~~4~ l;4j d~3 o0o6..... -5: -' 0.04 o -o - o-oe o.,o o.,z 0.,4, —- 0., 6 o.,.....o o.. ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^^Ei_\__ _____I^-_^^^NN^E~i.^^_ ^ I...... i 1-4~~~~~~~~~~~~~.~...r.....lt T~~~~~~~~~~~~~~~~~~~~~~~~~~ ^ o 5 --'^ —- --- --- -s>^^ |^,______^^i^^^^IIIW-^^i,,~~~~~~~~-"' fr~l.. ~: 1~~~~~~~~~~~~~........ Figure -6 ~. Temperature versus Time in a Un iform Slab with Unit Initial Temperature at All Stations...........;-t i:I d -i i If "S" t: i.. 0.5:~ W.............~~~~~~~~~~~~~~~~~~~~~~~~'l Ititi 1.. t...... I 1i;;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:: I: t~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~- r::: -iji~~~~~~~~~4 4 t t t f t 1 t + 0.5 1...,~~~~~~~~~~~~~..... imtd a,;::: irrrfrtr!;j~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~,t........... ~.. 0.4 T:: J, fIIr +.:~~~~~~~~~~~~~~~~~~~~~~~~~~ I I I U i: 1i I I j 1 I I ii i~ i 1f~ii i i t i 0.2........ i It 1f; t U I:: i i; f i:::,,: t~~~~~~~~~~~~~~~~~~~~~T fI it.......i: i I i, Tttii F fif f ti'f' I~~1 Ii~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~,f I i i.,..;.,,..1 t~~~~~~~~1 -F,, t~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Figur 2-6 Temeratue vesus ime i a Uifor Sla ~' Iwith nit IitialTempratur at Al Staions

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN for u2 (0) = 1 and the other temperature initially zero. In Figures 2-9 through 2-12 the same type of solutions are recorded for u3 (0) = 1, u4 (0) = 1, u5 (0) = 1, and u6 (0) = 1 respectively. If we denote the temper ature at the nth station resulting from a unit initial temperature at the mth station as u (t), then the complete solutions at the nth station for any nm initial distirubtion Um is simply N-1/2 u (t) -= m u (t) (2-12) m=l m nm For the case where N = 6-1/2, as in Figures 2-7 through 2-12, the temperature at the nth station is 6 u (t) = unm (t) (2-13) n m=l m nm Thus from the curves in these six figures the temperature distributions following any arbitrary intial distribution, approximated by 6-1/2 stations, can immediately be obtained. 27

I -I:- I-4Xi~~~~~~~i i k 1-l~ll.... 1-'' _ _ - _ _ _, ^ _ _ _ _ _ __ _ _ _ __.._._............. iiii~liiilii~il~ ti_.itiii-t — ii —i~f~f~- — 1 — ---- --- - -t-~' 1 -' 1' —-''1'- 1- - --- --- -'- r- -- - -T- - - - --- ---- - -1 *1 - ----- - 1:1-.-. - -- m - -. - t -- t - -- - _ _ — - - - -- - - -. - t t C -,.,t —,tt',, t'::tt;',:s','' —,;;'' 1:;-''-,1;1.......... — r ~-f o-'L... -. 1 -. -' -:' -. -.... - -': -'1-'.......................... t- l —-. -f. — - --— *. --- -..... < - 11: —. -— 1- 1 —-1 ---- ----- 1 ——: — 0 1...''`'` ~-,m, —, —-1-~,-~ ~~ * ~-I i- + A - -~-;-~~~+~ -..l.~-'....'.-. -'. I......... ~-~~ F I + —r — ~ -~ —t~-~~1'~': + < i 01 - 1 — 1 1 ------ 1 1 — i — I1+-I; i>tr- — li- ^...- -..|.. —........ + ji ti:..:I; -::.l-::::.::::l.:!+:: -... -..-:;::::.::.:.:.l.:.:1^... 4-tit'-'- + * t-'' +t- -— 1;- ~ 4.. -................................. - t -1......... i i::::i:: i i::i-:i,.... _:.:'::.....~:...............i.:........ i........' 1-:1 h+5-.;l cq. ++...''1'': - -''1'1-: - 1'': 1''1' — -- -— 1:' -' 1 -.j ~~- i'..+-.tE t:-~ — i-::1' 0'-'.-'.1,-~-.. —.....-..'+ S 1sX-' 1'- -~t~t. ~-~ —- ~~~~ ~ —-—. _.....T...s'''~~'~`~'~' ol,.. t,.- t- -.,........,.. _ _ j l. A.. - i i, T...............7.... It:t: t t I _ _ -' _ t-...... o o. o...................... O:1;..~~~~~~~~~~~~~~~~~~~~~~~~~~~~. ~~-~~~~~~~~~~~~~~~~~~~~~~~~~~~~- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ T':!~~~~~~~~~~~~~~~~~~~~~~~~~~~:...,,...,................,~ ~~~~~~~~ ~~ ~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~ ~-~~ ~-~-~~~~~-~~~~.~ -~-.1 —-~ ~~~1~~-~*~~~~~~~ —- ~-~- -~~ —- I 9 lim,, I-l-r; —ll lt-~~~~~~t — ~ C —, ~ — -~~~~~~~r-~-, —; 1~ -- ~-~....... -— ~~~~~~~~~ ——:.-......::::r:-n-:::1; ~~-~~~~- ~:;:T 2:;:1:2:!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -7~~~~~~~~~~~~" ~- ~.h ~ —-.~ ~~ rct —t+ —.. ~-tr-r -1.,-.-rr — -r —.C 0 1 -f1 — -- f fL ~ -~I -;1~I~~-1.-~:-n;::7: -~; 38m v83 dY 431 - $

CNJ.............................. - - -................................rl t~l~. Li~-~-!_ifi~-~ —-.~i-~iitit iiitriiiliiiii liil-iiilitiiiirijiiriii............ii~i iji~r~i~__tiIii~iliiiiliiii~iii~liiij~i.,::::::~:::~:::~::........:::............~::~::::::::: ~~~~~~~~~~....................... r w: ~ ~ ~ ~ ~ ~ ~~ 0 i~ ---- -::7:::......II1.....: tf I Ifilt~ o w _ —I —l;"'F_ -~- --— ~~~~~~~~~~~~~~~~~~~~~~~o,0 1:,:: 1 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~:: B I w~~ ~ ~ ~ ~ ~ a "~~~~~~~~~~'~.t~~~~~~~~.........,,~ ~~~-~- ~~~~~-~ —~-~- — ~t —- -- ~- ~...........~...... L c O U~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ =f~~~..... 0 03'~ I~ ~ t~ ~~~~ —~ I~ - t'~- c — t: t:;: "' d~~~~~~~~~~~~~~~~~~~~~~- -O ~ — -~e t fl~~~~~~~~~~~~~~~~~~~~~~~~E0.................~,,...,,., ~:; -'~~~~~" -- — ~~~~~~~~~~~~~~W-i~~~~~~~~~~~~~~~~~C ~-7 - - L - 7 -1-1 — -T1 -- +T - 0 0~~~~~~~~~~~~~ d 0Ua 38nilv83dlN31 z9- -~ tt- -1 c

. 1! e 1t t-ti!!~ - i ~l~L~~~~~~~~~~~~~~~~~~~~~~~~~~~-~i~~1~, I~: I~ +U:[:'-~.r Z - D I ~It T'' I''~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I1 t_'Tk-~~~~~~. Itl!... t. I I z -I. I....... +t,t+' +,,, +.,.,.e.,-..+. <-.. ___ _,.. _ —''' I 44-..-..,_.,.. ___ _ -— 11. Cd 11 1~:11-',,.. -.........- -...... -'' t.......''......,-:........'..............................I.......::' I:: I.....' I..........' I: I..._.. _ _.........''. 1_:''+"tt-rttrtct._~._,,.,-trt~ttr t'__f; < __ _.. -. - - - -'1'-1.'1:1. 1.N....,, i.............i._..:~- L~f~:0...... _' ~i'-'' —..''1:'''1'..'....:::::' I _ 1: ~~~~~~~....;:.... - - -- - - —................. 1 -....... -1.. -... l n l. 1.....~ ~~~~~~ ~~~~ ~~~ ~~~~~~~'-''....'-..':'''. -.. i -.-'."''.:1'1 1' 1'.. /1 11. 1~............ I I...... X-_-_~tte _t _I=f CTC- 1 —T 1~ r-f 1t +i- l 11 I:11t -..,... -..-. _-..,-............. 1 0:- -'.1:-::1: ~,.-_,........... -_........................t~~~tt~.~i~~- -.-., ~~~~~~~~~~~~~~~~te..i.-...... -.' _.'................X't1' <1 0'-! tii n...) r e D ln t n -................ o.... O —r- O:, ~. O. ~. ~ OOO........................ O.. I I.........................................., 30 / / d~~~~~~~~C' ~ ~ ~ ~ ~ ~ T —t-~.t~ I' —'.'1.'' -- r - -.:..:~ ~4C~': t. -I; CT: ""+'c''. ~;+`'r:T::::..,.,L/' _................. a a..............,... F q,.. ~~~~I' I~I I t It' I' I Il: iit::::::ii ~ ~I- ~~, ~ I I: I ~ __ I ~ I [ I i iI~..1 ~... ~~ LO~~~~~~~~~~~( d d O d~~~~~~~~~~~~~~~~~~~~C d c5 d5 o 38filV83dVN31 ~:II~i —~~ittrif~; - -~?- ~~~~;- ~~~- -3O-...

..................................... i~t.71. 1r j t~~~~i i,~ ~i ijii ii Ii fii::::;: t ttlli iiifiiI ~ii:1i~i ~ ~ _~l~~fi iii ~l ii fiiji i~i i o.. a-fflas ma- [~~~~~~~~~~~~~~~~~~~~~~::Iii::!!:;-w "T'^^^^^^^^^^^^^^^^^'^^""^^^^^^^^^^ —^-^-~~..... --- 77: ---. --- -—:: -- -- -- - - -'- -- -- ___ _ - J L __ t_ _; _ __ __ __ _ _ " * ":: t::' *:.:::: ": te l ^':.1 jft' -::: ^:::: ~ t~ I:: 1':::I:.:;'-:::::: 1: *::,::: -:: r: -::::r:::i::::::::I::::~:::I:::t::::t:::'::::'::::::::'::::':::::.':..................:................................................. i i i Ii:I II I:1 f f- -4.0 __ d_^ ^ nHn~ [n._^ _^ _ _]_i o d o o! o~ EEH ^^ EEE E:::^ 1^? ^ ^:I:.:::::':I::I::I:':.~:::.::::: - - - - - - - - - --- - - - - - - -- - -............... 3wniVd3dlM31 3!........... o E:.l1: lli tli::I:. I I I I.: I-: I i.................... 11 i. ~~~~~~~~~~~~~........ OD i I: E:~~~~~~~~~~~~~~~~~~~~~~:T;:I Ef~~~~~~~~~~~~~~~~~~~~~~~~:'::' i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~If1 I cde OD I: Ff r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~::1 1 1 la, a ~,:,, cc ~ ~ ~ ~ A 1 i:.r. 1 i o~~~~~~~~~~~~~~~~~P 1.:I hll 0 c) o (D LO d~ re C\11d 3 HnI V H13d M93

_....... ~- ~. --.... —. —;........ -.....t,-..t....t.....~ —:: _i~ t11 _! i —— ~ tii.1.E _il t i.z~ _:I _~1 ~ { 1 t l _i i [[ tt' 1~ eX - t keX.utu7,< V -E-7 I-S'Lf..:: t _ i. —......-.__t...t.. _._.....,...._.. _... —. - -- s....-...... _ _ _:::::~:: ~-; iil~j:- ~ ~-::'i:!I. -':Ctff=T.:~i'~ - "~ ~.i ~ ~.~; -~:!~ ~i!iiii ii liii;,;~:;~!ii~!1 iiiii! ~~~~~~~~~~~~~~~~~~~~~~~.. _'___. _-'<_._ _ _ __...... -,,:::.::-.t'i. tl~-.. _.., _, _..-',..t.+t_ - —.t _ X ~;T,;:_T r, ~i~:'~;-i~ ~.:!~~-:'i:~:~;~ ~-:.:.........-...-.: -:-::-.....:..._-:_ CIO~~;~ ~!! ~ ~ i ~!- ~ 4/1 ti,-L ~~~~~~~~~~~~~~~~~~~~0-I 0P-.......................~_. r I~~~~~~~~~~~fT ~ ~ ~ ~ ~ ~ ~ ~~~~~~P~-~~~~~~ E~' _ T < < < <~~~~~ t lt +t _ _ _ +t__,.nr ri 1 ++ rF- "f k... I:.... -...r |.....: t~~~~~~~t7 t a ~~~~i_i f i_ |. i i i. i.blS, + | i':,'. I! *E 4,, 111tjIfi, o a ~~~~~~~~~~~~~~::. T:.~:,i:::::,i ii,:::.i-i:::T::::.~ ~::~: ~:::::::::::::::::::::::::::::::~:::::::':::t'::::-:.': ~"::-~::~~ ~::~ -~::~:'::'.~ i::t:::il:::: iilt::::~~-::~ ~::~::::::::::::::::::::::::..... - -'-......::: ~''::::'"::................::......... "......... t o < X~~~~.......if fIi' - 0 e1.1- - --:'1 —1 1 0'.tX................................ o L _ _.._..... t............_.. _ _. _.. t!.:............... I.. 1 q::........:o~~~~~~~~~~~~.?:J...................'......!'i.................. 100001 —: {0 |00 |; 1-i:t |0-A —- llf 10000 0: E...........:-............... 1- ---- T t 7 -- - - - AH..............:2'::27::::::::'::i:',........:: - OD i:; I) C00 001Ae e ~ 77 - _7; [: 1:;:: 7'..:........::::~ o 32

......... - _- -. - --- ---- -. o...... l - - v- - s. F + - + _. _...._.._ ___.__. -t_...................................._. _ _ _._. CM -- 1-l- i~liilfii i:ilijii~i jiiiili iilli~ii~i'11~ii! i 1:f-iiiiiiiil! i~11i~~~:if il: ifillli:::i -!~ill::i ~ -1!~iii~if ii~ iiiiii~ 1~ ~11;li i t ~ i i:........ i.. _ -. t. - - -.... -: ITi I.... _ _ _.. t....__........... _ t... _~~~~~~~~~~~~~~~~.......... t-. —-— __.......==. t'' -'-t'j i..i.._f i~ if~i i f. isi i i i*i. t i I i r i iti i I _ i i — Iili <i ii t i i II I 1 j 1111 II I i i iii i. -.........OD-..........,.... -—, - 1 - _ -. -. i i - i i i i I. i I t I i i t If i i I i I ii 1 -~ 1 ~ i 1 - i t t;tl 4. ii i i - iil'i i Ui i iilf~ii............. -.::i - - - -I 1 — -.... - -:i.-. r i i -::I i i I i i i I i........................l. I t i. X 4-I -t;t,;;t.;. +..................t................fit:i.tIt: II A i~~~~t- i. 0-i —- t +:r -- w 4- iXl iliii iliii iiii iijl-:ti'ii 1 -' s - - * -''~~l~i~i ii-r- i('1(' _ -i~ Il | |1 11 EX[:1-Lt tit 1 -v t3X H t -- t- t f 1-f f 7 71X /U ETg L Ic.......................... ^U 4 ~~~~~~~~~~~~~~~~~~~~~r-I.::1~~~~~~~~~~I r,-'~' i:: I:: f::::t -.L!::':-: — i::! ~ ~', rf::!:]~.:-t-' L::::_:i! - ~..........:. —tl-........'_'.:: I::: t::: 1:+::......:-..:,.t-t t::.':-':: - _._:: -::...-:. ~::. —:1:_1:::.....:===.......... =.........=...= = = =......st Al -tw~t C tfe t~t L +.... —-w+-!-swt < —---— X r ff —|] l.......................: rm X X < 8+.. ——.;<.... C Y st r - i i t i 1 i.iIiI 1 t t r i i 4; i i.' i t i t. i.. i i i + w: I 1 —. 1 - - 1. T t:::: — +T: — t — t1 i i i tL 2 E' I i i < E I i ii;iii ii ii ii ii iii i.. jW.il I ~il i jli ii..'1,'~ ii;,"~,t't,':,,,- i:l iItl $ <t ] I. * — t- -.. t 4 -t t t _ t t't' r-,; t 4 i. -. e.... *... TiF U X I tt e W i~~~~~~~~~~~~~~~~~~~~~~~~~~lt W W rr~rt3tI t~_t4 SZD [<~_ i~ i. o 0).-. o ). -- -o:l~lnlv cl~j............. ~1+~~~~~~~~~~~~~~~~~~~~~~..........~.. r- ~~-.......... O O.......... — i-~t ~ i~ ~~ ~~- ~ ~ b i 4 -: ~ t OD rl - WOt t..; d -~ - ~ ~- d 33n 3d

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 3 HEAT FLOW IN TWO AND THREE DIMENSIONS; CARTESIAN COORDINATES 3.1 Heat Equation in Three Dimensions using Cartesian Coordinates From Equation (1-1) we have for the basic equation of heat flow in three dimensional Cartesian coordinates au 3 au a au C (x,y,z)- = - K(x,y,z)- + -K(x,y,z) — at 8x ax ay ay a oau (3-1) + -K(x,y,z) - + S(x,y,z,t) az az where x, y, and z are the spacial coordinates of the temperature u(x,y, z,t). Assume that we are concerned with heat flow within a rectangle of dimensions Lx, Ly, and L in the x, y, and z directions respectively. Then we define the dimensionless distance variables x, y, and z by x = (3-2) y y L y =a = a =- (3-3) L xy L xy L Y X Y z y L z - = a, axz (3-4) L L X z X z Following our procedure for the one-dimensional case in Section 1.1, we let the heat capacity C and the conductivity K be given by C (Lx x, Ly y, L z) Co C (x,yz) (3-5) and K (L x, Lyy, LZ Z) = K K (xz) (3-6) 34

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where OC and 0K are dimensionless and equal to unity for a homogeneous medium. As before we introduce a dimensionless time variable t given by K t = ~o t (1-10) C L2 o X From Equations (3-2) through (3-6) and Equation (1-10) the heat Equation (3-1) becomes au a au 2 C au 2 8 Ou C O-x- K-Sx- + + S(,y t) (3-7) OC at =ax OKax +xy ay K ay xz a xz OK where the heat source input S is given by 2 C L S (s,y,z,t) = - S(L x, L y, L z, - t) (3-8) K x y Z K K K 0 Boundary conditions on the temperature or the heat flux must be specified at the edge of the rectangular medium. To solve the equation by the difference method we must consider the temperature only at stations throughout the medium. Let Ulmn equal the temperature at the Ith station in x, the mth station in y, and the nth station in z. Further, let there be L, M, and N stations along the x,y, and z coordinates respectively. Then the heat balance equation at the Imn station is given by lCmnduImn ( UL IKi +l/2 m, n-Ulmn)-K 1n_(u - ul)j dt -~mn,' L' 1UM/Z,-m, n l11 + a2 M[2 [K / n(ul-, n - umn) - K1 m(,n mn - mln)] +a N2[KI, m,n+l/2(u,m,n+l - mn) - K m l n um.nl I., m, n-l/2 + Smn (t) (3-9) Boundary conditions are imposed in the same manner described in Section 2. 2 while initial temperature conditions must of course be specified at each station. The electronic differential analyzer circuit would in principle be similar to the circuit for one-dimensional heat flow shown in Figure 2-1. 3. 2 Solution by Separation of Variables in a HomogeneousRectangular Mediu As in the one-dimensional case we shall estimate the accuracy of the difference method by comparing the normal modes, shapes and decay 35

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN constants of the difference equations with those for the continuous medium. For an isotropic homogeneous rectangular medium with no heat sources Equation (3-7) reduces to au au au a2 2 2 - -it (3-10) - = 22 + ax z (3-10) at ax y X To solve this equation by separation of variables we assume that u (x,y,z,t) = X (x) Y (y) Z (z) T (t) (3-11) Substituting (3-11) into Equation (3-10), we have 1 dT(t) 1 d X() a 2 d ) a x 2Z(z) dXY ~ xz __a = _ ---- -= - 2 — - + 2J- --- *y- + 2 (3-12) T. dt X dx2 Y dy2 Z dz from which dT + PT = 0 (3-13) dt which has the solution T (t) = Ae- pt (3-14) Also from Equation (3-12) each of the terms on the right-hand side must be equal to constants -Px - ly, and -3z respectively, where P = + P+ P (3-15) * Px y z Thus d2X dX + x X 0 (3-16) dx2 d2Y P + 2 Y = 0 (3-17) dy2 a 2 Let us assume as boundary conditions that the temperature is everywhere and - + Z = 0 (3-18) dz a XZ zero on the boundaries. The solutions to Equations (3-16), (3-17), and (3-18) are then respectively 36

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Xi (x) = sin irx, i = 1,3,....(3-19) Yj (y) = sin jirx, j = 1,, 3,.... (3-20) and Zk () = sin krx, k = 1,2,3,.... (3-21) Hence the decay constant pik for the ijkth mode is given by 2 2 ( (2.2 2 23-22) ijk i +'xy) 3-2xz and the complete solution can be written as o00 o00 00 00 02 2 2.2 k2)t u(x,y, z,t) =Zl I A sin iwx sin jry sin krz e- (i +y 3 + xz )t i13 =l k= 1 ijk (3-23) The constants A. are evaluated by applying the initial condition 1jk u (x,y,z,0) = U (x,y,z) (3-24) and are given by A ijk = 8 f j U(x,y,z) sin irrx sin jrry sin kwz dxdydz (3-25) 0 0 For the case where the initial temperature distribution is a constant UO the constants Aijk become ijk 64U0 (3-26) Aijk j and ka odd (3-26) rr ijk = 0, all other i, j, k values 3. 3 Theoretical Accuracy of the Difference Method Let us now solve for the normal modes of the three-dimensional heat-flow problem when the difference method is used. From Equation (3-10) the difference equation at the Imn station becomes dumn 2 22 dt = L +llm.n-2Ulmn+Uul, mn)+ 2 MZ(ul m+1 -2mn+ut, m-l, + a Nz (um n+-2umn mnl) (3-27) xz 1 m n+l-mn I, m. n-1 where L, M, and N are the number of stations in the x, y, and z directions, respectively. As in the one-dimensional case, we assume the mode shapes 37

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN agree with those for the continuous temperature distributions. Thus let Ulmn = asin i]- sin sin kn e-Pijkt (3-28) Substituting Equation (3-28) in Equation (3-27) we find that P ijk= Z[L2(-+ M (1-cos L ) + 2 N2 (1cos ) M2(3-29) Expanding the cosine functions in power series we have i.2 2.2 2 2 2 2 i jir 2.2 k = 2 [ 2(1 2 + 2 )+ n i (I 2 +. 12L xy 12 12N (3-30) In the limit of infinitely many stations along each axis the ijk values given in Equation (3-30) clearly reduce to the exact values in Equation (3-22). Again we see that the error in the decay constant Bijkdecreases as the reciprocal of the number of stations squared. For the homogeneous rectangula medium considered above our assumption of exact agreement between mode shapes for the difference approximation and the exact solution has proved correct. 3.4 Solution in a Two-Dimensional Homogeneous Medium For two-dimensional heat flow in a square isotropic homogeneous medium with zero temperature on the boundaries it is easy to show that the complete solution is given by 00 ~~ 2 2.2 u(x,y, t) =_ 1j Aij sin iwx sin jrry e (i + )t (3-31) where Aij 4 J f U (x,y) sin irx sin jry dxdy (3-32) 0 0 For the case where the initial temperature distribution is a constant UO the A.. constants become ij 16U A.. = 2 ~, i and j both odd 3-33 A T/ij =, all other i and j values When the difference method is used the mode shapes are again sinusoidal, and the decay constant is easily shown to be 38

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN'ij = [LZ2(-cos ) + M2(1- cos (3-34) as compared with the value for the exact solution which is = 2(i2 +j2) (3-35) In Equation (3-34) L is the number of stations along the x coordinate and M is the number of stations along the y coordinate. The series expansion for Pij in (3-34) becomes 2 2j2 pi'= L'dl L" +... tj(-Jr ) (3-36) L I The difference equation is the same as Equation (3-27) with N = 0 and a = xy 1. Where the same number of stations are used along x and y coordinates, L = M and the equation at the Im station is dUim dT =u + u, - 4u +u + ut (3-37) dT 1= Ulm 1, m+l 4lm lm 1,m (m-13 where T L t (3-38) Consider the 7 x 7 cell problem shown in Figure 3-1. Let the walls be held at zero temperature and the initial temperature distribution be unity. Due to the symmetry of the problem we need only consider the equations for 6 stations. Thus from symmetry u21 = u12, u32 = 23 = u24 u43 = u33 = U34, u13 = u14. The six difference equations are dul u- = - 4u11+ 212 dT du12 - = Ul 1 412 + U13 +22 dT du13 dT - l 3u13 + 23 (3+3 dT (3-39) du 2 2U12 - 4U22 + 2u23 du23 d_. = U13 + U22- 3U23 + 33 39

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN du33 = 2u2 - 2u3 dT U23 - 233 The electronic differential analyzer circuit is shown in Figure 3-2. Recordings of the analyzer solutions for the temperature as a function of time at each of the stations are shown in Figure 3-3. / I - II - 2362 / — -— / I U11 U12 U13 U14 yU31 U32 U33 34 y X//// //S//-////- 1/'// ///.// U41 U42 U43 u44 / /// ////^/^^^.I I I Figure 3-1. Two-Dimensional Problem in Heat Flow. 40

TV *un-mpaiAI snoaua0ouuoH o-dojilosi u-e uT MOT, VeaH I-euolsuawuTa-oMjL ioj ^Tnoaio jaz6i-euV'2-7- ain2Tq 22ni2' —II- -1 91 1.^^^^ s'0 14 PI ^s| S'O

ttc, -+,.... ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~....,~ ct~~~~~~~~~~ -'T` -~. ~~-~ ~ ~ ~ ~ ~ ~ ~ ~........ F \ Q ^ ^ y i~ y ^ ^ t4...a.....t..............i:......:..^. r:::::: *:: ".:' ":-:; " **: *;:::;:::::::::::-'::;::-:;: — y^^ —-^ —^-~- ^ —--- ---— ~~ =' t — ^E —^^e- -^ - ^^ —H... -- —..,,~ctl~-r CCC.-i~ jt.. —- - HJ -c- -r --— ~ —^ C -? — S. -- -'i, - -- - -- - -- 4 ^05 — ^^ — ^1 ^ - — ^ -— ^ —~-1~ t~ —-(~-~~~~- i —^, -— ~ — -- --.- ft!5fttf c~c~ct~ect~crtht ~eet cC-t~te............ —. — - ^ ^ s -- ^s- -- Ltlt-t~.ttt~-Ctt-L~tt ec~~fi ~ i~fl~lCt~t —+-~- crtt —..........~t. ~1................ i_ -- - " " - -: -^ + —-:^l, —~ — - - -r s.- - — s - - - -s,- - - - - - - - - -- ^ - -- - - ^ — ^ ----- - - " ---? — - - ^ -- - s -- -^. ----- -.. -, -^ ^ ^^..... o n'r 1.::' I -- _ _ _ _ _ _ __~+~crt~i- ~.... — _ _...__ - -.. _ _l.~ ~t - + ^ - ^+ - ^-^.... F r 3-. A O. -t++t-rm —cct~ —+-~~~~~~~~~~~~~~~~~~~~........ + -~~~~~~~~~~~~~~~~~~ Cr. I - I..... 0.4~~~~~~4.............~ -C* —-~,.. C~. -I,.1.~.......... —.ffC-~- 1tf- t ----. t ( 0.2 0 I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~. 0. t. 0.2 0 3 0 4. 8 1~~~~~~~~~~~~~'" W~~~~~~~~~~~~~~~~~~~~~~~~~~~~, Fiue33 nlze ouinfrTo-iesoaIHa r)~~~~~~~~lw naSur stoicHooeeu eim

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 4 HEAT FLOW IN CYLINDERS 4.1 The Heat Equation in Cylindrical Coordinates In the previous chapter we considered the solution of the heat equation in three spacial dimensions using Cartesian coordinates, which are convenient when the boundaries of the three-dimensional medium are rectangular. When the boundaries are cylindrical in shape it is convenient to describe the spacial equation in cylindrical coordinates. Let us treat the case of a medium with constant conductivity and heat capacity throughout, although the electronic differential analyzer solution for heat flow in nonuniform media is in principle no more complicated. For a homogeneous isotropic medium the dynamic heat equation becomes from (1-1) au 2 C = KoV u + (4-1) 8t In cylindrical coordinates r, 0, and z the equation is given by atu- a u au- a 2u a + (rJ Zt) (4-2) C0 —". KI- 2F -2.- ^ ^ 1''+(.,'0 t) (4-2) 0 0 ~LFai ai r2 a0 aQz Assume the radius of our cylindrical medium is r0 and the height is z0. We shall then define a dimensionless radius r and height z by r = - - rr z = - = a -, a = - (4-4) z rz r rz z o0o o and a dimensionless time t by t (4-5) t 2t C r o at r ar r a rz z (46) 43

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where the heat source input S is given by r 2 C r2 S (r,0,z,t) = ~S (r0 r, O, z z, - t) (4-7) K K O O Boundary conditions on the temperature or the heat flux must be specified on the cylinder walls and the top and bottom. To solve the equation by the difference method we consider the temperature only at stations along r, 0, and z. Thus let utmn equal the temperature at the Ith station in r, the mth station in 0, and the nth station in z. Also let L, M, and N represent the number of stations along the r, O, and z coordinates respectively. Then the heat balance equation at the tmn station is given by du, mn 11+1 (mn = [+1/2,2 ul, mil, n (mn m-l, n dujmn L |(l+l/2)(u m numn) ( )(t-l/2)(utmn u m) (22r), ml, Zmn arz 2 i[u m n+lr 2lmn +l m n- lj mn(T) (4-8) where T = L2t (4-9) and q Imn(T) = S I Zm u T (4-10) ~i~mn* S,, -- (4-10) L L M N L Boundary conditions are imposed as described in Section 2. 2, while initial temperature conditions at each station must be specified. The electronic differential analyzer circuit for solving Equation (4-8) at each station, although somewhat involved, is still quite straightforward. 4. 2 Solution by Separation of Variable As in the case of Cartesian coordinates we will determine the accuracy of the difference method for cylindrical coordinates by comparing the normal-mode decay constants and shapes for the exact solution and the difference-method solution of the heat equation. For the exact solution by separation of variables we assume that u (r,S,z,t) = R(r)0(/) Z (z) (4-11) 44

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN which, when substituted into Equation (4-6) gives 1 dT 1 d dR 1 d2 a d'Z -p ==- r —- _ — + rz __ (4-12) T dt rR dr dr r 0 do Z dz" from which dT dt + PT = 0 (4-13) Equation (4-13) has the solution T (t) = Ae-t (4-14) where p is the usual decay constant. Since the first and second terms on the right side of Equation (4-12) are functions only of r and while the third is a function only of z, they must be equal respectively to constants -r3 -r0 and - z, where P= r+ Pz (4-15) Thus dZZ 13 + 2 (4-16) dz a rz and r d dR 2 1 d2 (— r — +r p (4-17) R dr dr r0 do The solution to Equation (4-16) is Zk (z) = B1 cos z +B sin - z (4-18) k rz ^ rz Here Pk are the discrete values which pz must take on for the boundary conditions on the ends of the cylinder at z - 0 and z = 1 to be satisfied, while Zk (z) are the corresponding normal modes in z. For example, in the case where the temperature is zero at both ends, only the sine solution in Equation (4-18) is admissable, and hence B1 = 0. Also in this case 4k/,rz = ki and 2? 2 Pk = arz k1, k = 1,2,3 (4-19) If, on the other hand, both ends of the cylinder are insulated so that the flux is zero at the ends, then only the cosine solution in Equation (4-18) applies and 45

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 3k = arz k 2r, k = 0,1,2,3,.... (4-20) Note that when the ends are insulated, pz may equal zero and a solution independent of z is possible. Returning to Equation (4-17), we see that both sides of the equation must be equal to a constant, say P 3, so that dZ0 d +p - = 0 (4-21) d02 P and ld r + ( ) -)R = 0 (4-21) r dr dr r The solution to Equation (4-20) is of course 0p (0) = D, cos4p. 0 + Dz sin P 0 (4-22) where pj are the discrete values which P3 must take on. Since 0(0) must be 3-9 rperiodic every 21r radians in 0, it is clear that PJ = j or 2. pj = j, j = 0,1,2,3,4.... (4-23) Substituting Equation (4-23) in (4-21) we have 1 d dR.2 )R= - -r-+ (r ) = 0 (4-24) r dr dr r If we let p = r. r, Equation (4-24) becomes l d dR.2 -d pdR + (l - ) R = 0 (4-25) pdp dp p which we recognize as Bessel' s equation of order j. We can therefore write the general solution to (4-24) as R (r) = E1 Jj (ipr r) + E2 Nj (4r r) (4-26) where Jj and Nj are jth order Bessel functions of the first and second kind respectively. Boundary conditions at the cylinder walls (r = 1) on the temperature or heat flux will determine the discrete eigenvalues which pr may take on. Since any one of these values will depend on the order i of the root... 46

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN of R (r) and the order j of the polar function P(0), we will denote the eigenvalues by pij. Thus from Equation (4-15) the discrete values ijk which the decay constant 3 may take on are simply ijk = ij + Pk (4-27) If we denote the normal modes in the r variable as R.., then the complete solution for an arbitrary initial temperature distribution can be written as 00 00 00 u(r z,t i=l j=0 k- Aijk Rij (r) j () Zk (z) e (ij+k)t (4-28) i=l j=o k=0 where the constants Aijk are obtained from the appropriate triple integral involving u (r, O, z, 0). For the axially-symmetric case, o =1 and j= 0 for j / 0. Thus the Bessel functions in Equation (4-26) are of Zeroth order. If the temperature is finite at r = 0, E2 in (4-26) must equal zero. Denoting the roots of Jo (p) as lX, X.... ki, we have for the complete, axially symmetric solution with zero temperature on the cylinder walls 00 00 u(r,z,t) Z Aik Jo () r)Z (z) e i +Pk)t (4-29) i=l k=O ik i k where the constants Aik are obtained from a double integral involving the initial temperature distribution u (r, z, 0). For the case where the temperature in the cylinder is independent of z, the solution with zero temperature on the walls is simply 00 u(r,t) = A J (X r) e i (4-30) i=l 1 0 1 where A. is given by10 A. f Jru(r,O J (X r) dr. (4-31) 1 Ll o2 i'0 4. 3 Theoretical Accuracy of the Difference Method for Cylindrical Heat Flow Now that we have written the exact solutions for the normal-mode shapes and decay constants for heat flow in a cylinder, let us compare these 47

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN with the values when the difference approximation is made. The latter are calculated by assuming that the temperature utmn at the Ith station in r, mth station in 0, and nth station in z can be written as Ulmn (T) = m Zn e T (4-32) i.e., that we can separate variables. Substituting ulmn given above into Equation (4-8) and letting 0ymn = 0 (no heat sources), we obtain the following equations: ^zNz Z+l - ( - 2 Pz) Zn Z-1 (4-33) rz l(2 —M ) m t+01 = 0 (4-34) and -( +/2) (R12 R1) - (1 - 1/2) ('h - Rl)j+(l l - ) RI - 0 (4-35) l - /' L- L 2 where P =Pr+Pz (4-36) The similarity with Equations (4-15), (4-16), (4-20), and (4-21) is evident. We already know from Section 2. 6 that the mode-shapes for Equations (4-33) and (4-34) agree exactly with those for the exact solution (i. e., are sinusoidal). If we recall that 0m must be periodic with period M, it follows from analogy with Equation(2-20)that the discrete values j which S must take on are given by 2 j.2 2 ] -M (1 - cos zj) = j2 Ll- J 2... T (4-37) 2TT M L 3M _ Note that in the limit of infinitely many cells M the %j given above reduces 2 to j, the value in Equation (4-23) for the exact solution. For zero temperature or heat flux required on the ends of the cylinder the discrete values Ok which z must take on are given by 2- kit 22k2 ki k 2N2 z 1 cs 2 k2 2 a( k i + ) (4-38) Nrz 1N 48

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Here again note that in the limit of infinitely many cells N the &k from above reduces to the value arz k ir given in Equations (4-19) and (4-20) for the rz exact solution. Having obtained in the case of the difference equations the formulas for the eigenvalues of the angular and height-dependent modes, and having demonstrated that the mode shapes agree exactly with those for the continuous equations, we are left only with the problem of comparing the eigenvalues and mode shapes for the radial-dependent modes. We know that the solutions Rij (r) for the continuous equations are Bessel functions of order j, but to complete the comparison we must solve the set of equations represented by (4-35). This is a formidable task, since it is easy to show that the mode shapes in this case do not agree with the exact Bessel-function solutions and hence it is necessary to recalculate and solve the characteristic equation for each number of stations L. Furthermore, the eigenvalue p. given in Equation (4-37) must be substituted in Equation (4-35), making the equations to be solved different for each j value as well as each L value. For this reason it was decided in this report to limit the consideration of the radial-dependent functions to the axially-symmetric case, i. e., the case for which j = 0. If j =0, p = 0 and Equation (4-35) becomes at the Ith station 1+7~ljiR = 0 (4-39) [(+ 1/2) (R+ - Rt) - (I - 1/2) (1~ - RL )+-Rt = 0 (4-39) This equation is the difference approximation to Bessel' s equation of order zero. Let us consider the case of finite temperature along the axis of the cylinder and zero temperature at the walls, so that the boundary conditions on Equation (4-39) require that R/2 - R_/2 = (4-40) and RL = 0 (4-41) Note that I takes on half integer values 1/2, 1-1/2, 2-1/2,.... because of the boundary condition at the axis. Thus Equation (4-39) is iterated L-l/2 times. The resulting set of L-l/2 simultaneous homogeneous equations can be solved for the L-l/2 49 _

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN eigenvalues 3io and the corresponding normal-mode shapes. The procedure for carrying out this computation is discussed in Appendix I. The eigenvalues io for the exact solution are equal to the square of the zeros i of Jo. Thus po,0 (2.4048)2 = 5.7831, p2 = 30.472, P3, = 74.887, P40 = 139.04, p50 = 222.93, P6, = 326.57, etc. The percentage deviation of the eigenvalues pio from these true values are plotted as a function of the number of stations or cells along r in Figure 4-1. For L greater than 10-1/2 cells the curves were extrapolated from the lower L values, which were computed to better than 0. 01 %. Note that the eigenvalue error falls off as 1/L2 as expected. In Figures 4-2 and 4-3 the mode shapes for a 4-1/2 cell solution by the difference method are compared with the first-four modes of the continuous J solution. Note the fairly tolerable agreement despite the crudeness of the 4-1/2 station approximation. In Appendix I the mode shapes for 10-1/2 cells are compared numerically with the Jo functions. 4. 4 Analyzer Solution of Axially-Symmetric Heat Flow with Initially Constant Temperature In the same manner we considered the one-dimensional heat flow problem in Section 2. 7, let us obtain here the electronic differential analyzer solution for the radially-dependent temperature starting with an initially constant temperature of unity. This corresponds to finding the temperature as a function of time in an infinitely long homogeneous and isotropic cylinder (or one of finite length with insulated ends) initially at unit temperature throughout and with walls held at zero temperature. The problem was set up on the electronic differential analyzer with 6-1/2 cells (L = 6-1/2) using an integrator time scale of 5 seconds. From Equation (4-8) we have at the Ith cell d = - [(I + 1/2) (u+ - ut) - ( - 1/) (u - u1)J (4-42) The boundary conditions of symmetry requires that ul/2 = U-1/2 while the boundary condition of zero temperature at r =1 requires that uL = 0. For 6-1/2 cells the equations become / = 2 (ull/2- -U/)] 50

IS'*OIdf 3eaH I-e3Tpu.llOj Jo sapoI/ ~uapuadaQ -IeTIpe oj o Sllao Jo jaqtunN jo uoTlounfi e se JuBesuoD X-oaaJ Jo uoTI'eTAO aeQ 9 uaOJca d *'I- aanjIT. S"1133 4o0 383wnN 001 08 09 OS 0t 0~ OZ 01 8 9 c ~ 2 I __ — -- I- O'001ETITT I I I - I --— 1- -- ------ 0'080'09/i / V / / /// / / /,,,//,,, o /I ~~~~~rfTVV ~f, /~-bC~0,oz / / - o'9- < H,!,-^^/ o -- o,_ —- ------ 80O'g- ~> - 90t0 -- ~3'0___ I_______!0 I'o-!Ii j~ ~~~~~~~~~~~~~~~,0 /,fii! i,~ ~.o~~~~~~~~~~~9,S'O 0I1','0/

I. 0.8 IST MODE 0.6 RI 0.4 0.2 --- CONTINUOUS SOLUTION 4-1/2 CELL SOLUTION 0.0 I.Ops. —---— I —--------------- I I I -- 0.0 0.2 0.4 0.6 0.8 1.0 r 1.0 0.8 0.6 0.4 R2 2ND MODE 0.2 0.0 -0.2 -0.4-0.6 I I I I 0.0 0.2 0.4 0.6 0.8 1.0 r - Figure 4-2. Comparison of 4-1/2-Cell Mode Shapes with Jo'Modes 1 and 2. 52

1.0 _0.8 _ -- CONTINUOUS SOLUTION 4-1/2 CELL SOLUTION 0.60.4R3 0.2 3RD MODE 0.0 -0.2 -0.4-0.6 I I I I 0.0 0.2 0.4 0.6 0.8 1.0 r o I.0 0.8 0.6 4TH MODE R4 0 0.2 -0.2 -0.40 -0.6 I I I 0.0 0.2 0.4 0.6 0.8 1.0 r Figure 4-3. Comparison of 4-1/2 Cell Mode Shapes with Jo' Modes 3 and 4. 53

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN -dU/ /- [2 (u2-l/2 - ul- /2) - (ul-l/ - U/z)] ^dU~i_ i 2_[3-1-/2 2 -1/)2 ( du= -1 [3 (u3_i/ u2_l/22) - 2 (u2 -l/ - ulil/,) (4-43) du3-1/2 2 / = _ [4 (u4 1/2 - u31/2 3 (u3_-/2 - 1/2)] du.l/2 24-1/2- = - [5 (S-1/ 2 -U4 _1/ 4 (U4-12 /2 -U dU5-1/2 2 r dTl/Z = +-L6 (U6-1/2 - U5-1/2) - 5 (U5-1/2 - U41/2) Analyzer solutions for the temperature at each station as a function of the dimensionless time variable t are shown in Figure 4-4. The exact solution for this case of unit initial temperature is given by10 u (r,) = 2 Jo(r) -2t u(r,t) = 2 e I (4-44) (=1 XIJl(\Xl) It is interesting to note that at t = 0 this series falls off considerably less rapidly than the equivalent series in Equation (1-34) for temperature decay in a slab. This in turn implies that the higher modescan play a somewhat more important role in considering heat-flow in a cylinder. In Figure 4-5 analyzer solutions for the temperatures are shown at each station for u1/2 (0) = 1, the temperatures at the other stations being initially zero. Figures 4-6 through 4-10 show similar plots, but in each case with the non-zero initial temperature being u1-1/2, u2 1/2' U3-1/2, U4 1/2, and u5 1/2 respectively. If we let u1i (t) denote the temperature at the Ith station when ui is initially unity, then the temperature u1 (t) at the Ith station for an arbitrary initial distribution U. at each station becomes L-l u! (t) = il u/ i (4t) U 54

0 0 l-l-. 1~L —--— ~~z~~.~.~_~~tmKn^~Jz~~iNz II-^ =- i 11r-1 1 - -I.. If -- --- --- i —- ------------— v- I-. -- -t^ -i —i.-l: I — -. -X I;- — V -K -i <1/I1Iici -.'I _i _',= i-t t'. t'-.,. I N 1:S/ = 0A I l o:.................0 ^ ====zr::===zrz=I.....:: 7......... 142WS;22WXHf HXXXX o m..:....:-. "4:': * ). * N * — * -- -..If: - - " "-.....*..... 1..-:i::: i I- I i - I - 1: - i-...........g.aa..~fflffiaa~^ffi~s^^inn^==4t..1.1.. 4.E~iiisiasffI^^ii^^^^zii~rz~dE-4:- CdI.7I: ~..:...;~.'7%::....:: -:::.';: l........:'':":'1;.. t...............:................................ ~ —-- p iK-y l-..^^ -- ^ —,,.? T....-::::,:;,.::::% i ~'........... 38 fnlV 3d N 31 55 4 tt -f 1 r-=3aniVU d 31 *7 rf::r: t I: rf r:: t ~~tr'-:'': i~''`5 5~''~'

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN For the curves in Figures (4-5) through (4-10) L = 6-1/2. From these six plots one can from Equation (4-45) obtain the time-dependent temperature at each station following an arbitrary initial distribution (or at least a 6-1/2-station approximation to the arbitrary distribution). 0 O o:: o^^ ^ oil o o o o oi —-— i::lll^lflllltll o':- - - - - - - - -_ -;_t................................ — _ _ _ j __ _ _ _::,'.:-~~ ~~.*:~~~~~~~~L..,'. I';''. "'; *' -..'t...... nn.... i.....I ii ig,,i......;:........ CM - - - - - - - -- --- -- -- -~.................. —......f................... i / i i i:'...:..:..:..^..:!........................... 0 0 o3 00 r,- (, if) ^ro o I, c- 0.........._______5 6 ___________'''~' ~ ~''''':I............ Ir:ll. lll'ii: ill~liii ll'li II: IIII:I' i I i i I i I " 0 i 1:: I i - -. I - I -: - I: 1 - I: - I - - i: I I. I I i I I I I I I:: I:Cd'' ~~~~~~~~~~''' I;: I~ ~ ~ ~ ~~~~~~~~~~~~~~~~~Q F'i I:i 0III Ia e I i/lllllll~~~l~~l -II I I I i I I: I II I i I~~~~~~llf 4 i liliilii::ii:11 III: I I I II IIII- I i.. i ~:I..~ i. IIII-:I- I'-I.. Q -:i:illili iif li ll -~1 lit llli lllll:: I I i i: I:: li:ili::l: il llli I I I!

............ 0.9 0^ I|r - - --- - - - — iliil ~iii ll:. - — i - - -- -liiiiilil - - -- -- -- - - - - - — I- ----- -................................................... LLJ ~~'Yr^^~~~~ —~ —-~~~~ —tT —--—._ 0.8 0.7 0.6 -~ - 0.5 - - V - - - - - - - - -- -- - - - -- - - - - - - - - - Qn ~ ~ Q iifi i n lii l lnl lfjfj i Iii E I:ni i - _^ _ _ _rlr z -zln _ _ _ _ _ _ _iizlz _ _ _ _ _ _ _ _ _ 0.4 0 i iiit lii i iiiifii ii iiiiji..........v.................... __: —~ __.....____ _.. _ ___ _ 0.0 0.02 0.04 06 Q08 0.10 0.12 0.14.0.16 0.18 0.20 a22 0.24 Figur e 4-6. Axially- Symmetric Tmrature istri bution with ~~~0.1~u /z o) ~:Il-iir i- -1100 - r: l1:t ll rllll II: ~ ll: r: ll lI i _i. ~ iii i1 iI I I: I: I III I II::i ItI ri. I i1 ~I:.I~l~ llrflll- l1:Illljll l~il~fj~ IIIIll:lllllj ~ I7 7 =7 u 1-1/2 (0) 1. ~ ~ ~ ~ ~ ~ 1.__L I _

1.0'- - if''.' -:'''!'........I...:.:iti 0.99:::i::-: __~__ ~ - __Il.:i(1]~I..........$... L-.... _ L _I..l —........ i.......:-......l-; —---.........:.........~........ 0.8'It'11i:i iii',''~ ~.. [C-:::!z.. i~ II:.:::,.:.........-..,.... ~'....':'': — I -t_ i I' ~ ~ ~'kI.'I I!: i ~: ": ii:L ~ j.'::,..,...... 0.7~ LlIJ::' ~~~ ~ ~ ~~'::'':''; "' ~~' L''': ~!~~~~:....:.:,.:",',:':':.. _U,,,...~..... ~~,,......:......~:~,,,:'ii — I s':~;...'.li:';::'. -- -- -i-: l/t /' —-"'- ii i~, I;I:: D~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0.5:: i ii!::~'. —.~~~~~~~~~~~~ —-I ——.- __. ""/ F t~ I ~_'~~ —' W ~~~~~~~~~~~~~~~~ —- -~- - w ~ ~~~~~~u~/ (o) -,t. J W 0.4 0.3...... - 11''~~~''''` ""'''''''~~~~' ~~ ~ ~ ~:I:.::I:: ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~ —~~~ —f::::...~~~~~~~~~~~~~~~~ 0.2 1...'''~ - - - f - 1 -I- 7 7 -''''''' 0.3t~~~~~~~~~~~~~~~~~~~~~~~~~~t-~~~~~~~~ ~ 0 0..4 006 00 0 0.2 01.1.8 02'' ~''~~'''~'''-""~''''' ~ Ii:'' IV2 --- C...] i i i.1 i''~ Fiue4 7 xal-Smerc eprtreDsrbuinwt'''''' ~ ~ ~ u(0

1.0 0.9 I 7::3....-:....:*:........:: *::: *':*"::*:: *::::. * *:. *. *:. I~~~~~~~~~~~~~<1:~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ w0.2,o Q.D 0.4IOZ 0.3 0.0 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 t s Figure 4-8. Axially-Symmetric Temperature Distribution with u3-1/2 (o) = 1.

1.0 i I _ _ 0.9 - - - r 4 0.8! { \ * -j i i 1 |: i _ A~ 0.71 i — k ~ -i::t _ C' a 0.65 \ t.4 — i- -~~-,- ---.. —-.- -- ---. —. I i I.:::: III\ -[._.-..L-__ -. --: __ ------- ---------- - -- --— I41 0.w I I i I {I;:~'-' — -, —---—,I 0.0 0.04 0.06 008 0.10 0.12 0.14 016 018 020 0 2 0 24 t _ Figure 4-9. Axially- Symmetric Temperature Distribution with I 4-1/2 —- ( -) = 0" W Fgr 4,9: I: S T: w i t h u4-1/2 (o) = ".

! iii^ 4 E ^':;: i ^::::;:::::!::if::;:::::::::::t............ ~T' Eiii:iilniIi i!ii fft Izz~ m._n~ m 0^~~ ~ ~~~~~~~~~~...., —'... f... \.,,;~,..I.........~;;...;:: t:'::-,::;:',;;:::::;::;'::::::!:.::~':I::.::"'::!.! —' ^:.::.:::::t:: _::: t "t:....'..........t.',-....' I I:'........~..............::~::-:? - _ _ ^ _ _ _ i _ _ _ _ _ _ __.........-.... z... t.._::with:.:.. (0):........:......:.:.I.....:.::::::::......!::::,:.:.::.:.....:..............~~~~~ ~ ~~~..!..!!,!:::!!::::::::::::~::::::!'! i!::! 1:::::::!:......:::::............~~~~~~~~~~~~ ~ ~ ~ ~................ 0. 6 Ii i ~ ~~~~~~~ ~~ ~~~~~....~! i!!!.:::........ i'......i~:...i i:....i: cr..... ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.....~~~~~~~~~~~~~~~~~ ~ ~ ~ ~ ~ ~ ~..j o,~~~~~~~~~~~~~W q....:. ON ~ ~ ~ ~ ~ ~ ~ ~ il:::':.: -::!..' i:.!. i. i.!~ ii i~'ii i:!:::".... -. 0.4~~~~~~~~ —- 0~~~~~~~~~~~~1.4:::.,::::::::2:::~i~:: i ~i:U 021 OD0 0.2 _ 0O4 0,6 0.8 0.10 0.12 0.14 O. 16 0.1is 0. 20 0.220,. Figure 4-10. A&xially-Symmetric Temperature Distribution with u 5_1/2. (0) =1.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 5 HEAT FLOW IN SPHERES 5.1 Heat Equation in Spherical Coordinates Next let us consider heat flow in spheres. As in the Cartesian and cylindrical cases we will investigate the accuracy of the difference method by comparing the normal-mode shapes and frequencies with those for the exact solutions. For a homogeneous and isotropic medium the heat equation in spherical coordinates r, 0, and 0 becomes from Equation (1-1) C au 1 a -2 au 1 a au 1 a2u - -- -(r ( — + — (sin 9 + + +S(r,,, t) K at r ar Or r sin 0 ae ae r sin 0 a3 (5-1) where r is the radius coordinates is the azimuth angle measured from the x axis in the xy plane, and 0 is the altitude angle measured from the z axis in the rz plane. Defining a dimensionless coordinate r by r r = -r (5-2) 0 where ro is the radius of the spherical medium and a dimensionless time variable t by t = K t (5-3) Cr 0o we have from equation (5-1) 1 a 2 8u 1 a 8u 1 8 u u = - -(r -)) + - (sin ) e ) + 2 S (r,, 9, t)(5-4) at r ar ar r sin 0 a8 ao r sin 0 8a where r 2 Cr2 S (r,0, 9,t) = ~ (r r, 0,, 0 t) (5-5) K 0 K The medium within the sphere can be divided into stations along r, 0, and S and the difference approximation to Equation (5-4) can be written in terms of the temperatures at these stations. Also Equation (5-4) can be 62

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN solved by separation of variables. The resulting equation in 0 is simply the one-dimensional wave equation [similar to Equation (4-16)] with sine and cosine solutions, the equation in 0 has as solutions the associated Legendre polynomials, and the equation in r is 1 d (r2 dR + ( + 1) 0 (5-6) r dr dr r where R is the radially-dependent temperature function, p is the decay constant, and I is an integer. We can compare the eigenvalues and normal-mode solutions from each of these equations with the equivalent eigenvalues and mode shapes when using the difference approximation. Since any temperature distribution within the sphere is a combination of the normal-mode solutions, we could therefore estimate the accuracy of the difference method as a function of the number of cells. Because the three-dimensional problem was treated in some detail for cylindrical coordinates in the previous chapter, let us consider here only the radially-symmetric case. From Equation (5-4) we have in this case au I a Z au -u = - 8(r -+ S (r,t) (5-7) at r ar ar Boundary conditions on the temperature or the flux must, of course be specified at r = 1. At r = 0 the symmetry condition requires that au/ar = 0. 5. 2 Difference Equation for Spherically-Symmetric Heat Flow Let us divide the coordinate r into N stations. From Equation (5-7) the equation of heat balance at the nth station is given by dun (n + 1/2)2 (u n+- un) (n-12 (un - unl)]+n ) (5-8) dT n L n where T = N t (5-9) and (n 1 0- (T) S (- 2 T) (5-10) n N NN The symmetry condition at r = 0 requires that Ul/2 - U-1/2 = 0 (5-11) - 63

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN The boundary condition at r = 1 will determine whether the total number of stations N is an integer or a half-integer. 5. 3 Solution by Separation of Variables in the Spherically-Symmetric Case Let us now solve the exact equation (5-7) with S = 0 (no heat sources) by separation of variables. The resulting decay constants and mode shapes will then be compared with the equivalent constants and mode shapes when the difference approximation is used. Assume that the solution u (r,t) is given by u (r,t) = R(r) ePt (5-12) Substituting (5-12) in Equation (5-7) and letting S (r,t) = 0, we have 1 d 2 dR r + PR = O (5-13) r dr dr If we let X (r) = rR (r), Equation (5-13) becomes d 2X d2 + P3X = 0 (5-14) dr which has the solution X = A cos P't + B sin Wt. Since R (0) must be finite, it is clear that A = 0 and that R (r) = Bsin; t (5-15) r For a boundary condition of zero temperature at r = it follows that = jir, or.2 2 Pj = 2 j = 1,2,3,.... (5-16) Thus the normal modes Rj (r) can be written as 1 Rj (r) = j- sin jrr, j = 1, 2,3,.... (5-17) where we have normalized Rj (r) so that Rj (0) = 1. The complete solution starting with an initial radial temperature distribution u (r, 0) = U (r) becomes 00 B. 2in j2r e2-18, u (r,t) -= sinjTrr e3 (5-18) ir j=l j where 64

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 1 B. = 2Tj/ r U (r) sin jrr dr (5-19) ~3 0 For the case where the initial temperature distribution is unity B. = -2 (-1)J. (5-20) J Note here that all modes are equally present at t = 0. 5. 4 Theoretical Accuracy of the Difference Method for Spherically-Symmetric Heat Flow For spherically symmetric heat flow in a uniform medium the difference equation at the nth cell is given in Equation (5-8). For comparison with our results of the previous sections, let 0n (T) = 0 and uN = 0 (i. e., zero temperature on the walls of the sphere). Letting un = Rn exp [-(j/N ) ] we have from Equation (5-8) at the nth station 1 + 2 1 2 1 2-21) (1 + R ) ) Rn + (1 n n-l Zn n+Zn Nn Rl= 0 (5-N1) where R-/ R1/2 and RN = 0 The resulting set of N-1/2 homogeneous algebraic equations can be solved for the N-l/2 characteristic roots pj and the accompanying mode shapes for each N. How this is accomplished is explained in detail in Appendix II, along with numerical results. The percentage deviation of the decay constants pj from the exact values given in Equation (5-16) is shown in Figure 5-1, where it has been plotted as a function of the number of stations N for the first four modes. Note again that the error falls off as 1/N according to expectation. For N > 10-1/2 the curves have been extrapolated from lower N values, which were computed to 0. 01%. In Figures 5-2 and 5-3 the cellular mode shapes are compared with the continuous ones for N = 4-1/2, while in Figures 5-4 and 5-5 the comparison is made for N = 10-1/2 stations. Comparison of these Figures with the corresponding ones for cylindrical coordinates indicates somewhat better decay-constant accuracy for a given N for spherical heat-flow but considerably poorer mode-shape agreement. 5. 5 Analyzer Solution of Spherically-Symmetric Heat Flow Starting with Unit Temperature As an example problem consider a sphere with unit initial temperature throughout and zero temperature at the walls. This problem was _________________________________ 65 ________________________

99'SllOD JO caqtLnN aoq. Jo uofounfi s SE MOIM BE3H OT3jauMnuS XII-eotaiqdS Joj Ju-3suoj-X3e3aQ jo uOTIeTIaa aOuaojadc'1-d ajnM- i S-1133 JO 38WBnN 001 08 09 09 Ov 0 1 8 9 ~_ _ _ I 0'0010'080'09-./-z-~- 0' 0~0-01- 0*9o'20'0z, O'-'< 9*09'0- -0/' I O'l- --,- - —..

1.0 0 0.8 IST MODE 0.6 RI 0.4 - CONTINUOUS SOLUTION 0.2- o 4-1/2 CELL SOLUTION 0.0 I I I I I I I I 0.0 0.2 0.4 0.6 0.8 1.0 r _ 1.0 0.8 0.6 0 0.4 R| 2ND MODE 0.2 0.0 -0.2- 0 0.0 0.2 0.4 0.6 0.8 1.0 r - Figure 5-2. Comparison of 4-1/2-Cell and Continuous First and Second-Mode Shapes for Spherically Symmetric Heat Flow. 67

1.0 - CONTINUOUS SOLUTION 0.8 08 \ o0 4-1/2 CELL SOLUTION 0.6 0.4 R3 3RDMODE 0.20.0 -0.2 - 0 -0.4 I I I I I I I I 0.0 0.2 0.4 0.6 0.8 1.0 r - 1.0 0.8 0.6 0.4 R4 4T\ MODE 0.2 0.0 -0.4 1, ~,,, 0.0 0.2 0.4 0.6 0.8 1.0 r -W Figure 5-3. Comparison of 4-1/2-Cell and Continuous Third and Fourth-Mode Shapes for SphericallySymmetric Heat Flow. 68

I.0 0.8 iST MODE 0.6 RI 0.4 CONTINUOUS SOLUTION 0.2 0 10-1/2 CELL SOLUTION 0.0 I, I I I i I i t s 0.0 0.2 0.4 0.6 0.8 1.0 r - lm 1.0 0.8 0.60.4 - R2 2ND MODE 0.2 0.0 -0.2 -0.0, 0.0 0.2 0.4 0.6 0.8 1.0 r o Figure 5-4. Comparison of 10-1/2-Cell and Continuous First and Second-Mode Shapes for Spherically-Symmetric Heat Flow. 69

1.0O --- CONTINUOUS SOLUTION 0.8- 0 10-1/2 CELL SOLUTION 0.6 0.4 R3 3RD MODE 0.2 0.0 -0.2-0.4 i 0.0 0.2 0.4 0.6 0.8 1.0 r 1.0 0.8 0.6 0 0.4 R4 \4TH MODE 0.2 0.0 -0.2 -0.4 0.0 0.2 0.4 0.6 0.8 1.0 r - Figure 5-5. Comparison of 10-1/2-Cell and Continuous Third and Fourth-Mode Shapes for Spherically- Symmetric Heat Flow. 70

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN set up on the electronic differential analyzer with N = 6-1/2 stations and a five-second integrator time scale. From Equation (5-8) the equations are du1/2 4 U2 -4/ -l/d = - 4 1-1(2/2 1- / 1J du-i/z 4[9 (5-22) d3-1/2 4 r 3- -= 16 (u41/2 - u3-2) - 9 (u3-1/2 - u2/ ) du_/ 4 4 dU4_1/2 = - 25 (u5-1/2 - U4-1/2) 16 (u4-1/2 - U3,1/2)] dT 81 L du5_ 4 /2 du5 /2 = — 36 (-u51/2)- 25 (u 5-/2 - u4/ /2) dT 121 Recordings of the time-dependent temperature at each of the stations are shown in Figure 5-6. 5. 6 An Alternative Method of Writing the Difference Equations for Spherically Symmetric Heat Flow An alternative method for writing Equation (5-7) for heat flow in a uniform sphere medium is au 1 a2(ru) au= (ru) + S (r,t) (5-23) at r 8r This suggests the following formulation for the difference equation at the nth station. du = -(n + 1) un+l 2nun (n -l)unl]t n (T) (5-24) ______________ 71

0.7-v ^- s-tt^ _ ~ _ ^4^4 — 1 ii - -:-^ I — 06- i —— l - ^i -I —— t — -I-~~- I —~t —~C~t~ — — ~ — 4 -- -~- - - ----- -— 4 —-- -- --- ---- - 6 — 0 I\ uCC-C-CI —-I-vi I -* v- V i'' 1Ki-0.4 0.2 __ Figure 5-6. Analyzer Solution of Spherically-Symmetric Heat Flow with Unit Initial Temperature and Walls at Zero Temperature.

- ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where T = N t as before. It is easy to show that the normal modes of Equation (5-24) with n (T) = 0 agree exactly with the modes given in Equation (5-17) for the continuous solution. The decay constants 3j are given by pj = 2N2(1 - cos J) (5-25) N which are exactly the same as the decay constants when the difference method is used for heat flow in a slab of unit thickness and both walls at zero temperature. 1314 Although the alternative expression in (5-24) for the difference equation exhibits no errors in mode shape, the decay constants in Equation (5-25) show considerably more deviation from the exact values than in our original difference formulation in Equation (5-8) (e. g., -0. 75% versus - 0. 23% for the first mode with 10-1/2 cells). Furthermore a representation similar to Equation (5-8) must be used when the conductivity K varies as a function of r, since in this case the heat equation becomes -au = —a F2 K(r)- +S(r,t). (5-26) ft r 3r 3a 73

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 6 CHANGE OF VARIABLE TO IMPROVE ACCURACY 6. 1 Regrouping of Station Locations Often a problem in heat transfer may be met where over the medium of interest the temperature changes very rapidly in one region and very slowly in another region. For small values of t the problem shown in Figure 2-4 is of this type. The temperature near the left wall decays quite rapidly at first, since the temperature gradients are very high in this region. On the other hand the temperature near the right wall decays imperceptibly at first. If, in solving the problem, we are particularly interested in the behavior near the left wall during the first few instants, then it is clear that the accuracy might be improved if we could space the x stations more closely in this region. This can be done without increasing the total number of stations N by letting a new distance variable y = x. If we rewrite the heat equation in terms of y and take equally spaced stations along y, this will in effect give us much closer spacing in the equivalent x stations. This is illustrated in Figure 6-1, which compares the station locations in the two cases when N = 6.5. 6.2 Solution of Temperature Distributions in a Uniform Slab The original equation to be solved is au aZu a- = u (1-14) at ax with boundary conditions u (0,t) = au(lt) 0 (6-1) ax and initial condition u (x, 0) 1 (6-2) We define the new variable y as y = x/2 (6-3) from which 74

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Ui U2 U3 U4 U5 U6 au =u U O2 U3 U4 UU U6 6 u8I0 L- 1-ax Figure 6-1. Comparison of Station Locations for the Space-Variable Transformation y = xx. a _a dy 1 -1/2 a 1 a (6-4) ax ay dx 2 ay 2y y Thus Equation (1-14) above becomes au = -. -a r i a^ (6-5) at - z y 2y -ay J with boundary conditions au(,t) (6-6) y (0, t) = a y (6-6) ay Next we divide the medium into N equally-spaced stations along y. The equation of heat balance at the nth station becomes from Equation (6-5) 75

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN du 1 N N N dun I r 1 (6 dT n Ltn n + - n) - n-. - (6-8)2 where a new time variauie T is defined as T = t (6-9) For the case under consideration, N = 6. 5 and the six simultaneous eyqations to be solved are du 1 u dT- 1 C - (Un i - u ) - -. —- f'' -2'-u ( U ( dT 2-. 5 l 5 du, 1 r 1 -- - = - -- (u ) U, - - y L (u- u3)- -3 dT 1.5 -2- 5 du 1 1 =.- L- (u - u,) 5 duI 1 dT 4 L4.5 5 4 35. 5 4 dT 5 L25.5 4.5 d6 _ -- -(u - _ dT 6'L 05 The initial conditions state that u1 (0) = u, u) = u3 (C) = U4 (0) = U5 (0) = 6 (0) = 1. This problem was set up or the elec.tronic differential analyze 76 andtherecrde teperturs a ea!,i E —ta,,.nare shwni Fiure4)

[tt i J 1 H Lt~t 3 t~f 4 M0W.tflt0-0t tl -I ttt|I t11-]*t i t}- Ikt 4.s W4Xll i t f~tCi ti ~ l~; ~~ ~ tt t 1 iLk r i-j: "'0';1; 0 o t~?~t~ ~ —- 1 i t 2 ii t i I, 14 4 i 1 J i -: OI IIIIII I IT T71T_7_ - 8 T T I Iz 1 I T T I T T T I.O~~~~~~. i6 t i X01<~~, Si12<.X fi1H U l WAL I I 1 T r 0.8 2 ttmt X mt <mX323SL3Li444titilljtt~lt Nit; ~F i- I t- a i I- t 1S iO mc + lift 4t t -g t X t t t4t tlS4ll1111 i+ Fttt~~ m? S~~t- uP.6 -4ti i O1 0.4 H+ I+ 0.2 I N11-~~~~~~~~~~~~~~~~~~ O 2 XX EX XXL <<Z I o=o F tm KtwEFH wttf~tt-'t~X~ilst~fl 0,04 0.2t- mI 1- 3 ttt tt3 t t[ < t * t f tt +t t+!i f|Figure 6-2. Analyzer Solution for Heat Flow in a Uniform Slab; Equal Stations along y, where y = /x.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN zero temperature. On the other hand Figure 2-6 shows the same solution when equally-spaced stations in x are used. In that figure u1 occurs at x = 0. 154 and the temperature decay is not nearly so rapid. In Figure 6-3 points from the computer solution in Figure 6-2 are compared with the exact theoretical temperature distributions obtained from Equation (1-34). Although the accuracy of individual points in the figure is not as good as the accuracies the difference-method displays in Figure 2-4, one can from the points in Figure 6-3 draw much more accurate temperature distributions for the first few instants of time t. This is due again to the closer spacing of the stations near the wall. There are other problems which have been solved by the author where a new independent variable not only grouped the stations at closer intervals where desired but also improved the decay-constant and mode-shape agreement with the continuous solution..________________ 78 _

1.0 0.000625 0.0025 0.8.0 -1 4 I / / 0.04 0.6 Q =0.10 SO~4 EXACT CONTINUOUS SOLUTION o POINTS FROM ANALYZER SOLUTION 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 DISTANCE x THROUGH THE SLAB Figure 6-3. Comparison of Analyzer and Theoretical Temperature Distributions Across Uniform Slab with Unit Initial Temperature.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 7 SOLUTION OF HEAT-FLOW PROBLEMS IN SEMI-INFINITE MEDIA 7. 1 Equation to be Solved Up to now we have only considered heat flow in media of finite dimensions. However, we may frequently be interested in calculating the temperature in media which are infinite or semi-infinite in extent. For example, consider the conducting bar shown in Figure 7-1. U=U' u-O u =O Figure 7-1. Semi-Infinite Bar. Let the left end of the bar be imbedded in a perfect conductor held at temperature U while the right end of the bar extends to infinity. Further, assume that the medium surrounding the bar is at zero temperature, and that the bar loses heat to the surrounding medium at a rate proportional to the difference between the temperature of the bar and the medium, i.e., proportional to u. This assumption, often called Newton' s law, is valid providing the difference between the temperatures, expressed in degrees absolute, is small. We will also assume that the thickness of the bar is small enough to allow us to neglect any temperature variations across the bar. Thus the temperature u within the bar is a function of distance x along the bar and time t, and the heat equation becomes 80

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN C (x) - K (x) -kV u (7-1) at ax ax Here C (x) is the heat capacity per unit length of the bar, K (x) is the conductivity, and kV is the energy lost to the surrounding medium per second per degree per unit length. Boundary conditions are u (O,t) = U (7-2) and lim u (x,t) = 0 (7-3) X -00 The second condition could equally well have required zero heat flux or temperature gradient at x = oo. Let us introduce a dimensionless distance variable x given by x = L (7-4) where L is a length yet to be defined. Also let us write C (x) and K (x) respectively as C (Lx) = C 09C (x) (1-7) and K (Lx) = Kc K (x) (1-8) where 9C (x) and OK (x) are dimensionless and where Co and Ko equal the maximum value of C and K respectively. From Equations (1-7), (1-8), and (7-4) Equation (7-1) becomes L2k c a8u a ( au v7 9(xIx Kx -, 8~K(x8X) K0 u (7-5). 2t) C L 0 If we let L k /K = 1, then v o ~~~L k=V^~ C(7-6) v Also define a dimensionless time t by K k t = o T = v t (7-7) CoL CO Finally, from (7-6) and (7-7) we can write Equation (7-5) as 81

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN C x au a au 7-8 at -x K a (7-8) For the case of a uniform bar, which we shall treat here, Oc (x) = K (x) 1 and au a 2u mu (xu-u (7-9) at 3x with boundary conditions u (0,t) = U (7-10) and lim u (x, t) = 0 X -+00 and initial condition u (x,0) = U (x). (7-11) 7. 2 Exact Theoretical Solution Before considering the solution of Equation (7-10) by difference techniques, let us examine an exact solution. For the case where the initial 15 temperature distribution is zero, one can show that the solution is given by I -x x u (x,t) = cosh x - ex erf (t + t ) +- e rf - t) (7-12) where we have let the boundary temperature UO = 1. In the steady state Equation (7-12) reduces to u (x,t) = e (7-13) which is clearly the correct solution to Equation (7-9) when au/at = 0. In Figure 7-2 the temperature distributions along the bar at various times t, as calculated from Equation (7-12), are plotted. Note that the temperature along the bar gradually rises to approach the exponential distribution given in Equation (7-13). 7. 3 Solution by the Difference Method Using a Change of Spacial Variable Since it would require an infinite number of cells of finite spacing to represent Equation (7-9) and the boundary conditions of Equation (7-10), it 82

0.0 0.8 cr0.6 i0.4 0.2 0.0 _ -____ 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 DISTANCE X ALONG THE BAR Figure 7-2. Temperature Distributions along a Semi-Infinite Bar at Various Times.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN is necessary to introduce a new independent variable y which will transform the infinite range in x into a finite range in y.' For example, let y = 1 - ex (7-14) In this transformation the boundary at x = 0 corresponds to a boundary at y = 0, while the boundary at x = oo corresponds to a boundary at y = 1. From Equation (7-14) a ~ dy _-x 8 a. _ d = eXa = (1 y)- (7-15) ax ay dx ay ay and Equation (7-9) becomes au a au u = (1 -y) (1- y)- - u (7-16) at ay ay with boundary conditions u (0,t) = 1, u (l,t) = 0 (7-17) -X Note that the steady-state solution e given in Equation (7-13) becomes simply 1 - y in the new variable y. Note also that the heat flux at y = 1 is zero even for finite temperatures, confirming the fact that our solution is independent of the boundary condition at x = oo (i. e., at y = 1). Now that the heat equation in (7-16) has only a finite range of spacial variable y, let us break the distance along y into N equally spaced stations. The difference equation at the nth station becomes dun i n1/2 n-/2 (1 - n) [(1 - +/ ) (un+l - n) - (1 - / ) (Un - un)- u dT N N n (7-18) Equation (7-18) is iterated N-1 times, where the boundary conditions require that u =, uN = 0. It is easy to show that the steady-state solution to the difference equations yields a temperature distribution which agrees exactly with the continuous steady-state solution e = 1 - y. The complete problem was set up on the electronic differential analyzer for N = 7 cells, and in Figure 7-3 the time-dependent temperatures at eaun station are shown. The analyzer solutions in the figure agree with the theoretical solution in Figure 7-2 to within the width of the recorder-lines. 84

-~~~~~~~......l.s........-.- - --...:..:......,.....-,. _ -.'.. -. _,.,,,,.l.-.,._. _._._..:1 -.;;; _.......................:.: +..; -..........::...:.,..,......::.................... r —tt-~ —r —~~~~~~~~~~~~~~ — -- -- _ ~~~~~~~~~~..........X-X x,__... _ _. _. _ -~~................. - H - -':1'.' 1... i.... * - - t.... tt-.......... - ~.. 4-i__ Li. _, *. _t-... —. _,+,.., w _ l_. _ w - i-__,__. [ —_e'I-...... ~-++ ~t- -t r-~+ = F G - - I t- - t- I t _ ttt — -r~.=tt t,. tIt-.-. >- - — t — -- t - - -t-~ — i- -- ~~~~~-~ —- I -T-+tt-_r t-t I - _ I I I m. W +> t _.. _._........ _...._, +_ - - -.....t t... _ - tt t_ I _t tt- -*-r ---. —*-+ — _ —mc$te t:-=.O; - -t:. t s: r- t~~~~~~~~~' i- tr i — ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ - *.'7.t' r"'. t T *;-+ = + +., T - - ei+ — -— 5^.-, - - + | _ _ 4- T 8 -+t _ _ _ - t..-., t-CL -14;t. _ _ ___f..___ +_|_ ___* t7 cn~~~~~~c;r'" t — * t.... — ~ —. -~ E - +__4-f........................... _ _ + _r — -_ EE X iE _ 4L- *-4 - = 1 -t=- 4t — Q _.*,. —-— t —t_-=_1@.......................Q)` + Sv-+ —T t -— r ___ --- t i + t - f - 1'. C —t- rl t-'.,. _.~-~ -t*C- 4 ltt 4 + tC t!t -_tt' _ + — -t t ___..__-t=t. =: — | -= —t f -- -— 1 —-+ _41ts' Iti @t-_ll 141'!+ 1|- _1 e _" + st+Q)-_e + _* __ _z|_++_ t|||* ___|_|I_:t T1! =e Cm *~~~~~~~~~~~~~~~~r..- + — rt t t __ - It + t _-~ =t r _ _ t tt - T tE — — h |r iSt= t l + - 0. t..twX.t__. XIIi\I, E+t i: *t'l t'i- * _cf _~~Z d - ili''"+''1'V tT_ t. — Ft'-:-.- -i-t. -t t-ftt- _ -' I' \ ~~~~~~~~?ctttm'-.,,, rl.l te —stc- - ~-..r-t- + t-.. V~- X1';' n c li o.-l-...t;~.' t -'-' ~~- t - *';._ _ =|= —_tt. -~ t t-'*-t —-- t_, I|,I'*t`'t-t'' t-'~tt —~t -~*-t -t + t _. - ir c H t W —7 —W-t — t~t= -~~ —l —-- ~ = — 1 —— ) i —— +ch-=Ii -|rI m Ot t_, _. _, r _ ~-. =t _ v ~-, t.:.. 0:t\ <=1,,1.' rrt.fs -; I-,: -1\ -Il'rt = 51l111l1,11 - -tt_, n,. = - 4 _. _:~~' LL.- I;i.n = i s I- t -tL, I,- _v|-,-,.W.=t — -- =,t —-, il_-_._,,;.IU m _l._.i.., - ( l l *,.i * * --- O- O O I~~,-_~m-~T —rctnti 3un udd ~ - -'~C~~ "t*t — 1 f+S~t ~-85 i t —

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN For the 7-cell example considered it is interesting to compare the equally-spaced station locations on the y scale with their equivalent locations on the x scale. This is done below. Equivalent Station y Coordinate x Coordinate 0 0. 000 0.000 1 0.143 0.154 2 0. 286 0. 336 3 0.429 0. 560 4 0.572 0.847 5 0.715 1.252 6 0.858 1.945 7 1.000 oo From the simple example considered in this chapter it seems clear that a wide variety of dynamic heat transfer problems through semiinfinite media can be solved by the difference method using the electronic differential analyzer. Without complicating the circuitry we could have included time-dependent boundary conditions, a non-uniform bar, and heat sources distributed arbitrarily through the bar. In each case the same spacial -X variable transformation, y = 1 - e, would be used. 86

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN CHAPTER 8 NONLINEAR HEAT-FLOW PROBLEMS 8. 1 Introduction All the heat-flow problems considered thus far in this report have been linear. For this reason, it was possible in most cases to treat the problems by separation of variables or the Laplace transform, even though the solutions were much more easily obtained with the electronic differential analyzer. However, in the present chapter we will treat a nonlinear heat-flow problem, one which will require the solution of a nonlinear partial differential equation. Solution of this type of problem by exact analytic techniques is not possible except in very special cases, so that the use of a computer of some type is virtually a necessity. The main difficulty in the theoretical treatment of nonlinear boundary-value problems is, of course, the failure of the superposition principle, i. e., the fact that the sum of two or more solutions is not itself a solution. 8. 2 Nonlinear Problem to be Solved As an example of a nonlinear heat-flow problem consider the temperature u (x,t) in the slab shown in Figure 1-1. Let the left side of the slab be hald at temperature Uo and the right side be insulated. Furthermore, let us assume that the conductivity K is a function of the temperature u, say proportional to u. Thus K = K u (8-1) and from Equation (1-2) the heat equation becomes C aU a [ K u a ] (8-2) oat ax 0 ax where for simplicity we have assumed uniform heat capacity CO and conductivity as well as no heat sources. We realize, of course, that nonuniform media containing arbitrary time-varying heat sources can be treated by the ___________________________________ 87 _________________________

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN electronic differential analyzer with essentially no added complexity. The boundary conditions are u(O,t) = uo u t) 0 (8-3) 0 BOx ax where L is the thickness of the slab. We denote the initial temperature distribution as u (x,O) = U (x). (8-4) Let us define a dimensionless temperature u equal to the ratio of the actual temperature divided by the initial temperature at the right wall. Thus let u (x,t) = u(xt) (8-5) U(L) Let us also introduce dimensionless distance and time variables x and t given respectively by x = (8-6) L and K U(L) Ko(L ) t = 0 (8-7) 2 C L 0 Equation (8-2) then becomes au 8 u at = ax-u au (8-8) af ax Uax with boundary conditions u (0,t) u, au(,t) = 0 (8-9) ax and initial condition u (x, 0) = U(x) (8-10) 8. 3 Formulation of the Difference Equations and Analyzer Circuit In the usual manner we next break the x dimension into N cells and write the equilibrium Equation (8-8) at the nth cell. Thus du dT = Un+l/2 (n+l - n) - Un-1/2 (Un -ni) (8-11) 88

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where T = N2t (8-12) We are immediately faced with the problem of representing the temperature u at half integral stations since the conductivity, which is a function of u, occurs at the half-integral stations. Let us approximate un+l / as the average of un+l and un Thus assume that u - u? n+l n \ unl/2 n+l n (8-13) from which Equation (8-11) becomes 2 2 2 du u ~ -2u + u dun Un+l - 2u + u n__" = _ ^ n ~n-~l (8-14) dT 2 The electronic differential analyzer circuit for solving this equation at the nth station is shown in Figure 8-1. Note that a servo multiplier, as described k kun_- 50 NOTE: k VOLT EQUALS + UNIT TEMPERATURE 100 I00 2 2 2 Figure 8-1. Analyzer Circuit for the Temperature at the nth Station; Conductivity Proportional to the Temperature. 89

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN in Section 1-5, is used to generate the square of the temperature at each station. In the Figure it is assumed that the temperature u never becomes negative. The time scale of integration is RC seconds. In the equation at the first station the prescribed boundary temperature Uo is introduced, while at the N-1/2 station we let uN+1/2 = UNj because of the zero heat-flow condition. N-1/z 8.4 Analyzer Solution for Uniform Initial Temperature Distribution Across the Slab The nonlinear heat-flow problem discussed in the previous section was set up on the electronic differential analyzer using 6-1/2 stations and with 50 volts equal to unit temperature (k = 50). Solutions were obtained for an initially flat temperature distribution of unity. Time dependent temperatures at each station are shown in Figures 8-2, 8-3, and 8-4 for leftwall temperatures U0 of 0, 0. 2, and 0. 5 respectively. Note that the temperature decays are not truly exponential, even after some time has elapsed. The actual type of decay which is present will be discussed in the next section. 8. 5 Exact Particular Solution by Separation of Variables; Comparison with the Difference-Equation Solution In order to establish the accuracy of our difference-approximation solution to the nonlinear heat-flow problem let us look for a particular solution to equation (8-8) by separation of variables. Thus let u (x,t) = X(x) T (t) (8-15) Substituting Equation (8-15) into (8-8), we have 1 dT _ 1 d dX = ( -X - -a (8-16) T dt X dx dx where a is a constant. Thus we have the two separated equations dT 2 dt+ aT = 0 (8-17) and d dX dXJX +aX = O (8-18) The time dependent solution T (t) is from Equation (8-17) 90

10 il^^^ST:~iitI11~1:rlf:: ~l:-ji -1 — --— _-F1 —----..........._ -------- - - - - ^*-^0.9:I~N^I^ —-Ir _ 0.8 I -f -— Is- n —I — 0.4 - 0.2 0.1 0.0 0,2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 8 TIME, t Figure 8-2. Analyzer Solution for Heat Flow with Conductivity Proportional to Temperature; Conducting Wall Held at Temperature U0 = 0. 91 O.S: -. 040..8j0. 14i.. TIMEP t ~ ~ ~ I-t~L Figre8-. nayzr oltin orHet lo wthCodutiit

i t ~ ~ ~ r -4ii ~;; - - ^ " 1^^ i| g^^^^^^^^^^^^^-^^p-^1 —----------— 0 ffi a) p~~~~~~~~~~~- "*i ^ -g|^^|gi||ll^-^ ^ ^ —^t —^ f -^-544~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~44- - TR:| |:,a;.t|| ^l | ^ ^ t ^ ^^ ^^f, t - 4-...._ 9 Z0. TTT~~~~~~~~~- -T T-j~~~~~~~~~~~~~~ ~...... 0. + +~~~~~~~~~~~~~~~~~~~ ----- ^'I iIli i4^fle i ^ - - -'^~ t T~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~~~~~~~t -1...... f...... Uj~~ 4 -4 o - 0~~~~~~~~r4 1 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ f T~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~f I i ^n I ^ I ^l~~iil~aili^^li-t|~ilill^-^i trt " Orr~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~r -~~~~~~~~~~~~~~~~~~~~~~~~J d- Tddd 4'aynivd1di~ia11 - 0 0 0 0 0 0 0 0 0 Q~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----- ~~~~i-l 1 ~ ~ ~ ~ f 8fi83bB

''El I ~tt st I, l iii~m^ ^^"^^~~~~~~~~~~~~~ I I...... i.i........ o.9w p J.:: 1 I: mj|| ||pttlIIItIIIIEIEIti t-tIIIIIIIIIII i 04 i. ^~F,......1...;^I.IIIlr0.37^1i Ijl~. *.....-'I " - - -— c —I — -' 02~~~~~~~~~~~~~~~~~~~~~~~~~.............. zm litf 0 1.-.... i ffs ni aiE=^sln: ==^=z^-=i4 00 111ii~iI!!111^1^^^1^ - L ~- 0-2 0-4 0.6 0.8|;01.2 14 — 16 18 TIME, t Figure 8-4. Analyzer Solution for Heat Flow with Conductivity Proportional to Temperature; Conducting Wall Held at Temperature U0 = 0.5. ~~~~~0.9 ~~~ ~~1 93 ~ ~ 44............. H't~~~~~~~~~t,~~~~~ii N t......... — t-~i~: 1~- I-:; I:::::II _ I: Ttt4, -~~- ~~L i I 1: ~ ~.N t 0.7,r 1: L i::t~~~~~t U i.......'' it r;i~~~ 1:....... 1J ijj 1;;....... tt I 1~~...... /ij i t~k~i i'~' 0.4 T t t i I~~~~~~~~.. p ti ~~~~~~~~~...... 0.3 i 7:1;i 0.2....trt,t; ~................. O i! i t-e-l i. i i-~~~~~~~~~~~~~~TIM, Fiue8-.AayzrSltinfrHetFo wt odutvt

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN T (t) = (8-19) at + 1 where we have let T (0) = 1. This is quite a different type of decay function than the exponential decays obtained in the linear case. For Equation (8-18) let us use the boundary conditions dX X (0) = (1) = 0; X(l) = 1. (8-20) This is a nonlinear eigenvalue problem which can be solved with the electronic differential analyzer. We let time on the analyzer be distance x through the slab, and the problem is started at x = 1. The integration proceeds toward x = 0. The value of a is varied for each solution until the condition X (0) = 0 is met. Thus several trial solutions are required to converge to the correct solution. The latter is plotted in Figure 8-5, and the value of a = 1. 114 was used for this solution. Thus if we start the problem with the temperature distribution shown in Figure 8-5 the distribution should maintain exactly this same shape while it decays according to Equation (8-19). To test our difference-approximation accuracy we set up the initial temperatures across the slab exactly as given in Figure 8-5. After release of the initial conditions the solutions shown in Figure 8-6 were obtained. Over the time interval shown the shape of the distribution stays constant to within one percent of the correct instantaneous values, and the time-dependent decay agrees with Equation (8-19) to within about one percent. For at least this special case, then, the difference-equation approximation gives a reasonably accurate solution to our nonlinear heat-flow problem for N = 6-1/2 stations. It seems hardly necessary to point out that we can solve problems where the conductivity is proportional to u, u, etc., by the same techniques used in this chapter. Problems where the conductivity or heat capacity are arbitrary functions of temperature can also be hafdled with the electronic differential analyzer by using function-generating equipment. 11 94

1.0 X(x) i / dx X dx + aX=O dX (i) 0.4 X( (O) d 0 X(l)= I a= 1.114 0.2 0.0 0.2 0.4 0.6 0.8 1.0 x Figure 8-5. Analyzer Solution X (x) for the Nonlinear Eigenvalue Problem. 95

[....,....::..... N.":"::" I e:.........................:,:':.: I:::':-' - f i.... I........ -....:."........ I.............. i.'..............:::'"......... 0........:...:......' —-..-.-.:::.:::: -::::::::::::.. 11 Iii {4 -I- -t- ~ ~'~~ I.'~'.... -\ft:: I 1 1:: 1 1 1 2 t 1 1 J, I - 9-6 - ~. [.. 0.7 -i......... -- - -- - - - --- - ---- - - - - - - --- --- 0.2 - ^ — -u - _ _ - -- -- _ _ _ _____ ____ 0.0 I__ ___'1~~2 34 TIME, t Solution by Separation of Variables.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX I CALCULATION OF DECAY CONSTANTS AND MODE SHAPES FOR THE DIFFERENCE APPROXIMATION TO AXIALLY-SYMMETRIC HEAT FLOW I. 1 Continuous Equations to be Solved In Section 4. 2 we saw that the radially-dependent temperature function R (r) in the case of axially-symmetric heat flow in a cylinder is given by the equation 1 d dR r Er d-r - + P R = 0 (I-1) If we are interested in the case of finite temperature at the axis and zero temperature at the wall, the boundary conditions become dR(o) R (1) = 0. (1-2) dr The solution to Equation (I-1) with the boundary conditions of Equation (1-2) is given in terms of Bessel functions of order zero. Thus the normal modes or eigenfunctions become 2 Ri (r) = J (ir), p io (1-3) Here the eigenvalues pio are given by Xi, where Xi is the ith root of J (p). i 0 I. 2 Difference-Equations to be Solved The difference approximation to Equation (I- 1) at the Ith cell is from Equation (4-39)'P io 1 (I + 1/2) (R+ - Re) - (I - 1/2) (Re - R_1)+- R 0 (I-4) where pio is the eigenvalue for the ith mode of the difference equation and 97

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN L is the number of stations or cells along r. The boundary conditions require that R 5 -R =0 5 (1-5) 0. 5 -0. 5 and RL = 0 (1-6) Equation (1-4) is iterated L-l/2 times. If we let 3io i 17 (I-7) the following set of L-1/2 simultaneous algebraic equations is obtained. 2R1. - (2 - i)R. = 0 Th1z.s-^ ^i~-^^1 = 0 1.5 2. 5 (2 -i) R1. 5+ 5 R 5 2.5 R3.5 (2 - i) R.5 + 2 5 = 0 i~1/2 R~ -( - iR-l + x/IZ R1=0~ (1-8) L 2 R -1 (2 -i) RL2+ L-2 -2 L-3 L -1.5 L- L-l L-2 I. 3 Solution by Means of Characteristic Equations The above set of equations will have a non-trivial solution only if the determinant of the coefficients vanishes. From this determinant we obtain a polynomial in Xi of order L-l/2, the roots of which are just the eigenvalues Xi. Thus for L = 1.5 the characteristic equation is 2 - Yi = ~ (1.5 stations), (1-9) --------- R 98 ---—.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN which has one root, y = 2. From Equation (1-7) the eigenvalue P!, o 2. 25Y1 = 4. 5, compared with the exact value of 5. 7831 for the continuous solution (see the next to last paragraph in Section 4.3). For L = 2. 5 we find for the characteristic equation 2 8 - 12yi + 3y = 0 (2. 5 stations) (I-10) which has two roots, Y1 = 0.84530 and Y2 = 3.1547. From Equation (1-7) p, 0 = 5. 2831 and 2,0 = 19.7169. These compare with the exact values for the continuous solution of 5. 7831 and 30. 472 respectively. Similarly for L = 3.5, 4.5, 5.5, and 6.5 we obtain 2 3 2 16 - 48 i + 30 yi - 5 Y = (2 - yi) (8 - 20 yi + 5 Yi ) = (I-1l) (3.5 stations) 2 3 4 128 - 640 yi + 720 iy - 280 i + 35 yi = 0 (4.5 stations), (I-12) 2 5 3 4 5 256 - 896 ~i + 3360 yi2 - 2240 Yi + 630 Yi - 63 Yi5 = (2 - yi) (128 - 896 yi + 1232 yi - 504 y3 + 63 i4) = 0 (5. 5 stations) (1-13) 438.85714 - 4608 yi + 11520 yi2 - 11520 Y3 + 5400 y4 - 1188 5 + 99 Yi = 0 (6. 5 stations) (I-14) In each case the roots yi can be found, from which the eigenvalues pio are calculated. The ith mode shape is obtained in each case from Equations (1-8) by substituting the known yi into the first equation, assuming R0 5 = 1, and solving for R1 5. The second equation can then be used to find R2 5' the third for R3 5, etc. 3. 4 I. 4 An Alternative Method for Obtaining the Mode Shapes and Frequencies Since the characteristic equations become very tedious to obtain and solve for more than 5-1/2 stations, an alternative procedure is used. Instead of solving for the characteristic equation we make a guess at the eigenvalue yi for L cells and use Equations (1-8) to solve for the R values at each successive station. If we have chosen yi correctly, RL should l_________________________________ 99 _

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN vanish. We can select R. 5 arbitrarily as equal to unity. Thus for 10. 5 stations the equations for the successive R values become from (1-8) R1.5 = (1 - yi/2) R0. = 1 - i/2 (R0. 5 R -= -0.5 + 0.75 (2 - yi) R15 R 5 = 0.666666667R1.5 + 0. 833333333 (2 - yi) R2 R 5= - 0.75 R2 5 + 0.875 (2 - yi) R3 5 5.5 2=-5 3-2 R = - 0.8 R3 5 + 0.9 (2 - yi) R4.5 (-15) 5.5 3.5 R6. 5= - 0.833333333 R4 5 + 0.916666667 (2 - yi) R5 R 75 = - 0.857142857 R5 5 + 0.928571429 (2 - yi) R6 R8 5 = - 0.875 R6 5 + 0.9375 (2 - yi) R7 5 R9.5 = - 0.888888889 R7 5 + 0.944444444 (2 - yi) R 5 R 105 - 0.90 R85 +0.95 (2 - i) R9 5 In general we will not have selected the correct Yi and RL (in this case R10. 5) will not vanish. However, if we have selected a yi close to the correct value (and this is easy to do by extrapolation from yi values for fewer cells), the RL value will be fairly close to zero and an extrapolation to a better estimate of Yi is readily made. After two or three trials the Yi value can be obtained to at least 0. 01% accuracy, and the ith mode shape has of course been calculated in the process of getting yi. The above technique was used to obtain the eigenvalues yi and hence pio for 10. 5 cells. Table I summarizes the percentage deviations of Pio from the exact decay-constant values as obtained from the squares of the roots Xi of J (p). A plot of the results and extrapolations from those results appears in Figure 4-1. In Table II the mode shapes for 4. 5 and 10. 5 cells are compared with the exact shapes given by Jo (Xi r). In each case the value of R. 5 was set equal to the theoretical J value at that point. 100

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN TABLE I PERCENTAGE DEVIATION OF DECAY CONSTANT WITH DIFFERENT NUMBERS OF CELLS FOR AXIALLY SYMMETRIC HEAT FLOW IN CYLINDERS Number Percentage Division of Cells 1st Mode 2nd Mode 3rd Mode 4th Mode 5th Mode 6th Mode 1.5 -22.19 2.5 - 8.65 -35.30 3.5 - 4.51 -19.60 -41.94 4.5 - 2.74 -12.28 -27.53 -45. 79 5.5 - 1.85 - 8.36 - 19.21 -33.06 -48.27 6.5 - 1.34 - 6.05 - 14.09 -24.72 -37. 03 -50.00 7.5 -19.09 8. 5 9.5 -19.03 10.5 - 0.51 - 2.35 - 5.59 - 9.98 -15.82 -22.42 11.5 -19.01 TABLE II COMPARISON OF 4.5 - CELL MODE SHAPES WITH Jo 1st Mode 2nd Mode Station Jo(k, r) R1 (J - R.) Jo(X2r) R (Jo R 0.5 0.9822 0.9822 0.0 0.9081 0. 908 0.0 1.5 0.8457 0.8458 -0.0001 0.3167 0.3088 0.0079 2.5 0.6012 0.6014 -0.0002 -0.2818 -0.2930 0.0112 3.5 0.2990 0.2993 -0.0003 -0.3621 -0.3739 0.0118 4.5 0.0000 0.0000 0.0 0.0000 0.0000 0.0 101

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN TABLE II (CONTINUED) 3rd Mode 4th Mode Station Jo( 3r) R (Jo - R ) J(X 4r) R (J - R ) 0.5 0. 7819 0. 7819 0.0 0. 6147 0. 6147 0.0 1.5 -0.2185 -0.2658 0.0473 -0.4008 -0.5293 0.1285 2.5 -0.2381 -0.2554 0.0173 0.2675 0.3764 -0.1089 3.5 0.2879 0.3219 -0.340 -0.1304 -0.1873 0.0569 4.5 0.0000 0.0000 0.0 0.0000 0.0000 0.0 TABLE III COMPARISON OF 10.5 - CELL MODE SHAPES WITH J 1st Mode 2nd Mode Station J(Xr) R1 (JR) RJ ((- R 2r) RI (Jo - Rl) 0.5 0.9967 0.9967 0.0 0.9828 0.9828 0.0 1.5 0.9707 0.9707 0.0000 0.8505 0.8502 0.0003 2.5 0.9197 0.9197 0.0000 0.6126 0.6118 0.0008 3.5 0.8457 0.8457 0.0000 0.3167 0.3153 0.0014 4.5 0.7516 0.7516 0.0000 0.0204 0.0184 0.0020 5.5 0.6410 0.6410 0.0000 -0.2211 -0.2235 0.0024 6.5 0.5181 0.5182 -0.0001 -0.3673 -0.3699 0.0026 7.5 0.3878 0.3878 0.0000 -0.4003 -0.4027 0.0024 8.5 0.2548 0.2547 0.0001 -0.3276 -0.3294 0.0018 9.5 0.1239 0.1239 0.0000 -0.1794 -0.1804 0.0010 10.5 0.0000 0.0000 0.0 -0.0000 0.0000 0.0 3rd Mode 4th Mode Station 0(k3r(J- R) J) Rl (Jo R1) J(x) - RI) 0.5 0.9580 0.9580 0.0 0.9227 0.9227 0.0 1.5 0.6529 0.6508 0.0021 0.4069 0.3990 0.0079 2.5 0.1892 0.1842 0.0050 -0.1881 -0.2012 0.0131 3.5 0.2185-0.2253 0.0068 -0.4008 -0.4121 0.0113 4.5 0.3997-0.4060 0. 0063 -0.1599 -0.1617 0.0018 102

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 3rd Mode 4th Mode Station J (x r) R3 (J - R) Jo( 4r) R1 (J - R) 5.5 -0.3128 -0.3163 0.0035 0.1962 0.2035 -0.0073 6.5 -0.0560 -0.0556 -0.0004 0.2883 0.2965 -0.0082 7.5 0.1973 0.2009 -0.0036 0.0631 0.0642 -0.0011 8.5 0.3001 0.3046 -0.0045 -0.2011-0.2071 0.0060 9.5 0.2094 0.2123 -0.0029 -0.2202-0.2267 0.0065 10.5 0.0000 0.0000 0.0 0.0000 0.0 0.0 5th Mode Station 0J(X 5r) RI (Jo - RI) 0.5 0. 8776 0.8776 0.0 1.5 0. 1479 0.1307 0.0128 2.5 -0.3871-0.4097 0.0226 3.5 -0.1851-0.1885 0.0034 4.5 0.2431 0.2582 -0.0151 5.5 0.2111 0.2199 -0.0088 6.5 -0.1459-0.1552 0.0093 7.5 -0.2206-0.2313 0.0107 8.5 0.0668 0.0713 -0.0045 9.5 0.2146 0.2257 -8.0101 10.5 0.0000 0.0000 0.0 103

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX II CALCULATION OF DECAY CONSTANTS AND MODE SHAPES FOR THE DIFFERENCE APPROXIMATION TO SPHERICALLY-SYMMETRIC HEAT FLOW II. 1 Continuous Equations to be Solved For spherically symmetric heat flow the radially-dependent temperature function R(r) is, according to Equation (5-i3), the solution of the equation 1 d 2 dR - (r2 dR + R = 0 (II-i) r dr dr For the case of finite temperature at the center and zero temperature at the wall the boundary conditions become d = R(1) = 0 (II-2) From Section 5. 3 the exact normal-mode solutions Rj to Equation (II-i) are given by R(r) = -L sinfjr (11-3) Rj (r) jrr where the eigenvalues pj representing decay constants are given by.2 2 Pj = J2 = r,, 3,.... (11-4) II. 2 Difference Equations to be Solved The difference approximation to Equation (II-1) at the nth station is from Equation (5-2) 1 2j 1 2 (1~ + R (2 + (1 ) Rn+l - + — ) Rn + (1 - ) Rn- (II-5) 2n 2n N 2n where R. 5 R05 = 0 and R = (II-6) ________ 104

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Here N is the total number of cells or stations into which the radial coordinate is divided, while Pj is the eigenvalue for the jth mode of the difference equation. Equation (11-5) is iterated N-l/2 times, giving the following N-l/2 simultaneous algebraic equations: 4R1.5 -(4 - Yj) R0.5 0 16 R 2 ( - yj) 5+R0 = ~0 9 2. 5.j 1.5 0. 5 36 R 52 16 0 25 3. 5 - ( - Yj) R2.5 2-5 R15 (II- 7) n- + 0. 52 8n2 + 2 n-05 2 ( n 5 ) Rn ( 2 - j) Rn + (- ) Rn = 0 n n1 4n n n 2 2 2N- RN - N 2 N-i 5 L (N2) +2 i.. = 0 ( ) NR N2J N- 21N-3 N-2 N-i 4(N- 2)2 N-2 *- *( - ) 8(N- l)2 - RN + (N - 1.5RN-2: 0 N 4(N - ) - N - 1 where - j (11-8) Yj = N(II-8) 3j N2 II. 3 Solution of the Difference Equations As explained in Section I. 3 the above set of simultaneous equations will have a non-trivial solution only if the determinant of the coefficients vanishes. Solution of this determinant will give us the characteristic equation, an N-l/2 order polynomial, the roots of which are the eigenvalues y. Once the eigenvalues are found the mode shapes are obtained from the following equations, derived directly from (11-7). R1 5 = (1 - 0.25 yj)R R2 5 = - 0.25 R0 5 + (1.25 - 0.5625 j) R 5 (II-9) i05

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN R3.5 - 0.444444 R + (1.44444 - 0.694444 yj) R R5 = - 0.5625 R2 + (1.5625 - 0. 765625 yj) R3 5 etc. Here we can let R0 5 arbitrarily be equal to the correct value for the continuous solution, as obtained from Equation (11-3). Since the determinant of the coefficients in Equation (11-7) is very tedious to solve for large values of N, it is easier to use the technique discussed in Section I. 4 to evaluate the eigenvalues yj. In this method we guess at yj for N cells, substitute into Equations (11-9) and find RN. If RN does not vanish, we make a second calculation with an improved yj, and in this manner quickly converge to the exact yj values. This method also gives us the mode shapes and was used to prepare Table III, which presents as a function of N the percentage deviation in decay constants p j from the exact decay-constant values given in Equation (11-4). The results are plotted in Figure 5-1. Table IV compares the mode shapes for 4-1/2 and 10-1/2 cells with the exact mode shapes given by Equation (11-3). The results are shown graphically in Figures 5-2 through 5-5. TABLE III PERCENTAGE DEVIATION OF DECAY CONSTANT FOR DIFFERENT NUMBERS OF CELLS FOR SPHERICALLY-SYMMETRIC HEAT FLOW Number of Cells 1st Mode 2nd Mode 3rd Mode 4th Mode 5th Mode 1.5 -8.90 2.5 -4.47 -25.38 3.5 -2. 19 -16.28 -33.58 4.5 -1.29 - 9.75 -25.48 -38.01 5. 5 6.5 -0.61 - 4.60 -12.34 -23. 19 7.5 8. 5 9.5 10.5 -0.2 13 -.73 4.73 - 9.11 -14.74 106

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN TABLE IV COMPARISON OF 4.5 - CELL MODE SHAPES WITH R. = sin jwrr jTrr 1st Mode Station R (- r) Rnoe R1R R 0.5 0.9798 0.9798 0. 0 0.9207 0.9207 0. 0 1.5 0.8270 0.8619 -0.0349 0.4i35 0.5157 -0.1022 2. 5 0. 5642 0. 5992 -0.0350 -0. 0980 -0.0960 -0.0020 3.5 315 0.3135 0.2823 0.0312 -0.2015 -0.2506 0.0491 4.5 0.0000 0.0000 0.0 0.0000 0.0000 0.0 3rd Mode 4th Mode Station R3(r) Rn R -R R4(r) R R2 Rn 0.5 0.82699 0.82699 0.0 0.70532 0.70532 0.0 1.5 0.00000 0. 1517-0.151i7 -0.20675 0.14700-0.05975 2.5 -0.16540-0.29574-0.13034 0.09207 0.03962 0.05245 3.5 +0.11814 0.17697-0.05883 0.03499 -0.01042-0.02437 4.5 0.00000 0.00000 0.0 0.00000 0.00000 0.0 TABLE V COMPARISON OF i0.5 - CELL MODE SHAPES WITH Rj = sin j7rr 1st Mode 2nd Mode Station R (r) R R - n Rr) Rn Ri Rn 0.5 0.99623 0.99623 0.0 0.98516 0.98516 0. 1.5 0.96676 0.97398-0.0072 0.87103 0.89849-0.0275 2.5 0.90932 0.91948-0.0102 0.66658 0.69901 -0.0324 3.5 0.82699 0.83822-0.0112 0.41350 0.43953-0.0260 4.5 0.72410 0.73519-0.0111 0.16113 0.17516-0.0140 5.5 0.60598 0. 61605-0. 0101 -0.04528 -0.04394 -0.0013 6.5 0.47865 0.48708-0.0084 -0.17487-0.18309 0.0082 7.5 0.34841 0.35480-0.0064 -0.21723 -0.22990 0.0127 8.5 0.22150 0.22568-0.0042 -0.18301 -0.19463 0.0116 9.5 0.10370 0.10567-0.0042 -0.09909 -0.10569 0.0066 10.5 0.00000 0.00000 0.0 0.00000 0.00000 0.0 107 _____

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 3rd Mode 4th Mode Station R3(r) Rn R - R R4(r) R R- Rn 0.5 0. 96676 0. 96676 0.0 0. 94139 0.94139 0.0 1.5 0.72410 0.78129 -0.0572 0.54308 0.63497 -0.0919 2.5 0.34841 0.39767 -0.0493 0.04981 0.09333 -0.0435 3.5 0.00000 0.01525 -0.0153 -0.20675 -0.23135 0.0246 4.5 -0.19356 -0.20882 0.0153 -0.14517 -0.18336 0.0382 5.5 -0.19748 -0.22242 0.0249 +0.04478 +0.04072 0.0040 6.5 -0.07437 -0.08844 0.0141 +0.12819 +0.15179 -0.0236 7.5 +0.06445 +0.06851 -0.0141 +0.04834 +0.06298 -0.0146 8.5 +0. 12778 0.14247 -0.9147 -0.06686 -0.07708 0.0102 9.5 0.09169 0.10339 -0.0117 -0.08187 -0.09822 0.0164 10.5 0.00000 0.000000 0.0 0.00000 0. 00000 0.0 108

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN BIBLIOGRAPHY 1. Howe, C. E., Howe, R. M., Application of Difference Techniques to the Lateral Vibration of Beams Using the Electronic Differential Analyzer, E.R.I. Report 21i5-1-T, OOR Contract No. DA-20-018ORD-21811, Feb. 19.54, 2. Howe, R. M., Electronic Differential Analyzer Solution of Beams with Nonlinear Damping, E. R.I. Report 2115-2-T, OOR Contract No. DA-20-018-ORD-21811, April, 1954, 3. Howe, R. M., Theory and Operating Instructions for the Air Comp Mod 4 Electronic Differential Analyzer, E. R. I. Report AIR-4, ONR Contract N6 onr 232 (23), April 1953. 4. Hagelbarger, D. W., Howe, C. E., and Howe, R. M., Investigation of the Utility of an Electronic Analog Computer in Engineering Problems, Report UMM-28 of Aeronautical Research Center, April 1949, 5. Howe, C. E., Howe, R. M., and Rauch, L. L., Application of the Electronic Differential Analyzer to the Oscillation of Beams, Including Shear and Rotary Inertia, Report UMM-67 of Aeronautical Research Center, Jan. 1951, 6. Corcos, G. M., Howe, R. M., Rauch, L. L., and Sellars, J. R., Application of the Electronic Differential Analyzer to Eigenvalue Problems, Project Cyclone Symposium II on Simulation and Computing Techniques, May 1952, 7. Howe, C. E., Howe, R. M., and Rauch, L. L., Solution of Linear Differential Equations with Time-Varying Coefficients by the Electronic Differential Analyzer, Project Cyclone Symposium II on Simulation and Computing Techniques, May 1952. 8. Howe, R. M., Propagation of Underwater Sound in a Bilinear Velocity Gradient, Engineering Research Institute Report AIR-3, ONR Contract N6 onr 232 (23), March 1953. 109

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 9. Howe, C. E., Howe, R. M., Solution of Linear Differential Equations with Variable Coefficients by the Electronic Differential Analyzer, Transactions of the Institute of Radio Engineers, Professional Group on Electronic Computers, vol EC-2,(1953), 10. Churchill, R. V., Fourier Series and Boundary Value Problems, McGraw-Hill. (i941). 11. Korn, G. A., Korn, T. M., Electronic Analog Computers, McGrawHill ( i952). 12. Salvadori, M. G., Baron, M. L., Numerical Methods in Engineering, Prentice -Hall. ( 1952). 13. Haneman, V. S., Solution of Partial Differential Equations by Difference Methods Using the Electronic Differential Analyzer, E. R. I. Report AIR-1, Oct. 1951. 14. Haneman, V. S., Howe, R. M., Solution of Partial Differential Equations by Difference Techniques Using the Electronic Differential Analyzer, Proceedings Institute Radio Engineers, vol 41 (1953). 15. Churchill, R. V., Modern Operational Mathematics in Engineering McGraw-Hill (1944). _____ 110 _____.