ENGINEERING RESEARCH INSTITUTE DEPARTMENT OF AERONAUTICAL ENGINEERING 2115-4-F UNIVERSITY OF MICHIGAN ANN ARBOR SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS BY DIFFERENCE METHODS USING THE ELECTRONIC DIFFERENTIAL ANALYZER Office of Ordnance R1e-search, U. S. Army Contract No. DA- Z 0-1 8 -ORD- 2181 1 Final Report, Period from December 1, 1952 to May 1, 1954 Submitted for the Project by R. M. Oewe Associate Professor Department of Aeronautical Engineering

SUMMARY This final report summarizes the investigation of the solution of both linear and nonlinear partial differential equations by difference techniques using the electronic differential analyzer. Complete details of the research are given in three technical reports. Two categories of physical problems are considered:, (1) the lateral vibration of beams, and (2) heat flow. In both cases time is preserved as a continuous variable, while the spacial variable is broken into stations or cells. Thus spacial derivatives are approximated by finite differences, and the resulting set of simultaneous ordinary differential equations are solved by the electronic differential analyzer. Theoretical accuracy of the difference method as a function of the number of cells along the spacial variable is investigated by comparing the eigenvalues (normal-mode frequencies for the beam, decay constants for heat flow) and normal-mode shapes with those for a continuous medium. Results show suprisingly few cells are required for representation of the first few modes to several percent accuracy. Theoretical as well as computer solutions are obtained for cantilever, hinged-hinged, free-free, and clamped-clamped beams, both uniform and non-uniform. Beams with viscous -damping and time-varying boundary conditions were also solved, and the vibration of beams including deflection due to transverse shear are treated by the difference method. Cantilever beams with nonlinear damping terms such as velocity-squared damping and Coulomb damping are solved successfully. Dynamic heat-flow problems treated by the difference method, both theoretically and with electronic differential analyzer solutions, include one, two, and three-dimensional flow in rectangular media as well as flow in cylindrical and spherical media. Change of independent variable to improve accuracy and to solve flow in semi-infinite media is demonstrated. The nonlinear heat equation for a medium having a conductivity

proportional to temperature is solved with the electronic differential analyzer and results are compared with a particular exact solution. As a result of these investigations it seems clear that the electronic differential analyzer is an excellent tool for rapidly solving many partial differential equations, both linear and nonlinear. There is a one-to-one correspondence between resistor values in the circuit and physical parameters in the problem, and the desired dependent variables (e.g., beam displacement, velocity, bending moment, or temperature, heat flux, etc. ) are all available as time-varying output voltages of the computer.

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. (AD ATTN: ORDTX-P Department of the Navy Washington 25, D. C. Office, Chief of Ordnance Washington 25, D. C. Commander ATTN: ORDTX (Technical Library) U. S. Naval Proving Ground Dahlgren, Virginia Chief, Detroit Ordnance District Director 574 East Woodbridge Director Detroit 31, Michigan Applied Physics Laboratory Detroit 31, Michigan Johns Hopkins University Commanding General 8621 Georgia Avenue Aberdeen Proving Ground, Md. Silver Spring 19, Maryland ATTN: BRL Corona Laboratories Commanding Officer National Bureau of Standards Frankford Arsenal Corona, California Bridesburg Station Philadelphia 37, Penna. U. S. Naval Ordnance Lab. Philadelphia 37, Penna. White Oak, Silver Spring 19, Commanding Officer Maryland Picatinny Arsenal ATTN: Library Division Dover, New Jersey The Director Commanding Officer Naval Research Laboratory Redstone Arsenal Washington 25, D. C. Huntsville, Alabama ATTN: Code 2021

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 ATTNim RDR eMa dOffice of the Chief Signal Officer ATTN: RDR 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. 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

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.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS BY DIFFERENCE METHODS USING THE ELECTRONIC DIFFERENTIAL ANALYZER INTRODUCTION Whenever we wish to solve partial differential equations by means of the electronic differential analyzer, it is necessary first to convert the equations to ordinary differential equations, since the analyzer can integrate with respect to only one variable,> namely time. If the partial differential equation is linear, this conversion to ordinary differential equations can often be done by separation of variables, which results in ordinary differential equations of the eigenvalue type. The normal modes, or eigenfunctions, can then be found, after which the complete solution is built up by combining the normal modes. The above method of separating variables and obtaining a series type of solution can be carried out fairly efficiently on an electronic differen.1-4 tial analyzer. Certainly, for most problems the analyzer is much faster than any hand methods. But for the engineer who is interested in getting quantitative answers to specific problems even the analyzer approach might seem somewhat tedious. Then too, the control engineer would often like to have a real-time simulation of the system being controlled, and in many cases partial differential equations are required to describe this system. It therefore would be highly advantageous to solve the partial differential equations directly with the electronic differential analyzer. This can be done by replacing some of the partial derivatives by finite differences in order to convert the original partial differential equations into a system of ordinary differential equations. Assume we are interested in solving a partial differential equation in which the dependent variable y(x, t) is a function of both a distance variable x and a time variable t. Instead of measuring the variable y at all distances x, let us measure y only at certain stations along x; thus, let Yl -----------------. —--------------

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN be the value of y at the first x station, y2 be the value of y at the second x station, yn be the value of y at the nth x station. Further, let the distance between stations be a constant Ax. Now clearly a good approximation to y/ax l / 2 (i. e., the partial derivative of y with respect to x at the 1/2 station) is given by the difference 3y Y, - Yo (11) (1-1) ax 1/2 Ax In fact the limit of equation (1-1) as Ax4 0 is just the definition of the partial derivative at that point. Writing (1-1) in more general terms ay Yn - Yn-1 ax I1/2 Ax In the same way 8y (ax x n/ (1-3) / ax I Ax ax n^+1/2 ax n-/I or from equation (2-2) y ________(1-4) 8x (AYx)2. (1-4) n Thus we have converted partial derivatives with respect to x into algebraic differences. The only differentiation needed now is with respect to the time variable t, so that we are left with a system of ordinary differential equations involving dependent variables yo(t) Y, t),... y SOLUTION OF THE HEAT EQUATION 2. 1 Basic Equations for Heat Flow The basic equation of heat flow is given by au c 6 = V.KV u +S (2-l) at where u = temperature and is a function of the spacial coordinates and time, 2 __________________________

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 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 (2-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 (-KV u) is a vector representing the heat flux. The components of -KV u 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 -KV u, as well as the initial temperature distribution throughout the medium. 5 In technical report AIR-10 the form of Equation (2-1) and its solution by difference methods and separation of variables is discussed for Cartesian coordinates, cylindrical coordinates, and spherical coordinates. For example, in Cartesian coordinates for heat flow through a homogeneous slab we have the equation 3u. a2u au = 8u (2-2) at ax where x is dimensionless distance through the slab and t is dimensionless time. For boundary conditions u(0, t) = uo(t) (2-3) (lt) = 0 (2-4) corresponding to a prescribed temperature uo(t) at the left boundary and zero heat flow at the right boundary, the difference equations become 3 —------------------- 3

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN dU1 1 du 1 r [U u l + uot) dt (Ax) dt u3- 2u2 + ul dt (Ax)2 N-l/Z 2UN-3/2+UN-/zj.. (2-5) dUN-1/2 1 dt(-)2 I UN-1/2+N-3/22 lT+UI "^ —' (Ax U N- 1 / 2 + UN-.2 dt (Ax) L N1I2N3I2, where Ax is the interval between stations along x and N-1/2 is the total number of x stations across the slab. This set of N-1/2 simultaneous first-order differential equations can be set up and solved directly on the electronic differential analyzer. In the technical report5 the time-dependent solutions for the temperature at each station are shown for a number of representative initial conditions. 3. SOLUTION OF THE BEAM EQUATION 3. 1 Basic Equations for Lateral Beam Vibrations The basic equation for lateral displacement y of a beam is given by 2 2 92 a2 EI() - + () = f(x,t) (3-1) a t2 OR ax at where x = distance along the beam, EI = flexual rigidity, in general u function of x, p = mass per unit length, in general a function of x, t - time f(x, t) = external applied force per unit length. The bending moment M is given by 4

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 2 ay (3-2) M(x,t) = EI(x) - (3-2) 8x while the shear force V is OMn2 (k) t) V(x t) a M(,t) (3-3) Boundary conditions depend on the type of end fastening. For a cantilever (clamped end at x = L y(L,t) = ay(L,t) = 0. (3-4) 8 x For a free end at x. L M(L,t) = V(L,t) = 0 (3-5) while for a simple supported (hinged) end at x = L y(Lt) = M(L,t) = 0. (3-6) In technical report: AIR-7 the normal-mode frequencies and shapes for Equation (3-1) are given for uniform free-free, cantilever, hinged-hinged, and clamped-clamped beams, along with theoretical solutions using the difference approximation. Electronic differential analyzer circuits and solutions using the difference technique are presented for the above cases and for nonuniform beams, beams with vidcous damping, beams with time varying boundary conditions, and beams for which Equation (3-1) is modified to include transverse shear effects (this is necessary for beams whose thickness is not small compared with their length). In addition, technical report AIR-8 presents analyzer solutions for beams with nonlinear damping terms, both velocity-squared damping and coulomb (dry-friction) damping. 3. 2 Difference Equations for Uniform Cantilever Beam As a simple example, consider the equation for lateral displacement y of a.uniform..beam. It can be written as 4 2 + yA - f(xt) (3-7) ax4 at -------------- 5,

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where x is dimensionless distance (the beam length in x is unity) and t is dimensionless time. To solve Equation (3-7) using the difference method we consider the transverse displacement y only at equally spaced stations approximately as mn = | (Yn + 2Yn + n - 1 (3-8) where Ax is the distance between stations. In the same way the complete difference equation at the n-th station becomes d2 d -2~ n(m2 =+ -2m +mn + n +(t) (3-9) dt (Ax.)n n-1 where fn(t) is the applied force at the n-th station. For a built-in end at the 1/2 station, the boundary conditions of Equation (3-4) imply that YO = YI = ~ (3-10) while for afree end at the N + 1/2 station the boundary conditions of Equation (3-6) imply that mN = N+1 = 0. (3-11) Thus the complete set of difference equations for an N-cell cantilever beam becomes d Y1 2 2 d Y3 2 d2y3 1 2- = - ( 4 - 2m3 + m2) + f3(t) dt (Ax)2

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN d YN12 1 N - - 2= - -2 (mN - 2m + mN 3)fN-(t) dt (AX)2 N- dt2 ( 2 - + mN -2) + fN - (t) dyN 1 d —t = - - ( mNl) + fN(t) dt (2 x) N1 2 where 1m = 2 (Y2) (Ax) 2 m2 = - (y3 - 2y2) "3 Z (Y Y2) m3 = 2 (Y-4 - 23 + Y2) (Ax) mN-2 2 (YN -1 2N- 2 + YN - 3) (Ax) 1 mN-1 -( YN - 2N-1+ YN-2) Thus we have converted the original partial differential equation given by (3-7) to a set of N - 1 simultaneous second-order ordinary differential equations which can readily be solved with the electronic differential analyzer. Computer output voltages represent 7

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN time-dependent displacement, velocity, and bending-moment at each station, while the external forces at each station are represented by timevarying voltage inputs. For computer circuits and recorded solutions of this and many other representative beam problems, the reader is referred to the technical reports. 6'7 For convenience the detailed list of topics and figures covered in these reports is given in Appendices II and III. 8

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX I Outline for Technical Report AIR-7 Entitled APPLICATION OF DIFFERENCE TECHNIQUES TO THE LATERAL VIBRATION OF BEAMS USING THE ELECTRONIC DIFFERENTIAL ANALYZER TABLE OF CONTENTS Chapter Title Page PREFACE i ILLUSTRATIONS v 1 INTRODUCTION 1 1.1 Basic Equations for a Thin Beam 2 1. 2 Solution by Separation of Variables 5 1. 3 Replacement of Partial Derivatives by Finite Differences 8 2 DERIVATION OF THE DIFFERENCE EQUATION FOR BEAMS 11 2. 1 Equation for the n-th Cell 11 2. 2 Boundary Conditions 13 2.3 Initial Conditions 15 2. 4 Summary of Variable Transformations 15 3 ELECTRONIC DIFFERENTIAL ANALYZER CIRCUIT FOR SOLVING THE BEAM EQUATION BY THE DIFFERENCE METHOD 17 3. 1 Linear Operations of the Electronic Differential Analyzer 17 3. 2 Analyzer Circuit for a Cantilever Beam 19 3. 3 Analyzer Circuit for a Hinged-Hinged Beam 19 4 THEORETICAL ACCURACY OF THE DIFFERENCE TECHNIQUE FOR UNIFORM BEAMS 24 4. 1 Uniform Hinged-Hinged Beam 24 4. 2 Uniform Cantilever Beam 26 4. 3 Uniform Free-Free Beam 31 4. 4 Summary of Theoretical Accuracy Calculations of the Difference Techniques for Beams 34 9 - --

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Chapter Title Page 5 APPLICATION TO CANTILEVER BEAMS 35 5. 1 Transient Response of a Uniform Cantilever Beam 35 5. 2 Static Deflection for a Uniform Load 35 5. 3 Determination of Normal Mode Frequencies 37 5.4 Component Accuracy Requirements 40 5.5 Effect of Voltage Transients 41 5.6 Tapered Cantilever Beam 42 5. 7 Uniform Cantilever Beam with Concentrated Mass Load at the Free End 46 6 APPLICATION TO HINGED-HINGED BEAMS 51 6. 1 Analyzer Circuit for the Hinged-Hinged Beam 51 6. 2 Uniform Hinged-Hinged Beam 52 6. 3 Uniform Hinged-Hinged Beam with Concentrated Mass at the Center 53 7 APPLICATION TO FREE-FREE BEAMS 56 7.1 Uniform Free-Free Beam 56 7. 2 Stability of the Free-Free Beam Circuit 58.8 BEAMS WITH VISCOUS DAMPING 61 8. 1 Beam Equations Including Viscous Damping 61 8. 2 Difference Equations Including Viscous Damping 62 8. 3 Computer Circuit for the Viscous-Damping Case 63 8.4 Impulse Response of an 8-Cell Cantilever Beam with Viscous Damping 64 9 BEAMS WITH TIME-VARYING BOUNDARY CONDITIONS 68 9. 1 Method of Introducing Time-Varying Boundary Condtions 68 9. 2 Cantilever Beam with Specified Displacement at the Free End 68 9. 3 Beam on Elastic Foundations 69 10. VIBRATION OF BEAMS INCLUDING DEFLECTION DUE TO TRANSVERSE SHEAR 71 10. 1 Equations for Transverse Beam Motion, Including Shear Displacements 71 10. 2 Difference Equations Including the Transverse Shear Effects 73 APPENDIX I - CALCULATION OF MODE FREQUENCIES AND SHAPES FOR CELLULAR FREE-FREE BEAMS A-i APPENDIX II - CALCULATION OF MODE FREQUENCIES FOR CELLULAR CANTILEVER BEAMS A- 12 10

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Page APPENDIX III - MODE FREQUENCIES AND SHAPES FOR CELLULAR CLAMPED-CLAMPED BEAMS A-15 BIBLIOGRAPHY A-1 5 ILLUSTRATIONS Figure No. Title Page 1-1 Cantilever Beam 3 1-2 Station Arrangement for Cantilever Beam 9 3-1 Operational Amplifiers 18 3-2 Analyzer Circuit for a Cantilever Beam 20 3-3a Response of 8-cell Uniform Cantilever Beam to a Uniform Impulse; Stations 6, 7, and.8 21 3-3b Response of 8-cell Uniform Cantilever Beam to a.Unifonm Impulse; Stations 2, 3, 4, and 5 21-a 3-4 Hinged-Hinged Beam 22 3-5 Analyzer Circuit for Hinged-Hinged Beam 23 4-1 Normal-Mode Frequency Deviation for a Uniform Hinged-Hinged Beam 27 4-2 Normal- Mode Frequency Deviation for a Uniform Cantilever Beam 29 4-3 Comparison of Mode Shapes for Cellular and Continuous Uniform Cantilever Beams 30 4-4 Normal-Mode Frequency Deviation for a Uniform Free-Free Beam 32 4-5 Comparison of Mode Shapes for Cellular'and Continuous Free-Free Beams 33 5-1 Bending-Moment Response of 8-c1 Cantilever Beam to a Uniform Impulse 36 5-2 Excitation of the Second-Mode of an 8-cell Cantilever Beam 39 5-3 Tapered Cantilever Beam with Concentrated Mass Load 43 5-4 Displacement at Wing-Tip Following a One-Second Uniform Impulse 44 5-5 Comparison of the Moments at Stations 1 and 4 of the Tapered Cantilever Beam with and without Concentrated Mass Load 45 5-6 Uniform Cantilever Beam with Concentrated Mass Load at the Free End 46 11

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Figure No. Title Page 5-7 1/3 L versus ML/MB for an 8-cell End-Loaded Cantilever Beam 50 6-1 Uniform Hinged-Hinged Beam with Concentrated Mass Load at the Center 53 7-1 Cellular Free-Free Beam 56 7-2 Drift of an 8-Cell Free-Free Beam Following Release of Zero Initial Conditions 59 8-1 Analyzer Circuit at the nth Cell Including Viscous Damping 64 8-2 Unit Impulse Response of a Uniform Cantilever Beam with Viscous Damping (Input and Displacement at Stations 7 and 8) 65 8-3 Unit Impulse Response of a Uniform Cantilever Beam with Viscous Damping (Displacement at Stations 2, 3, 4, 5, and 6) 66 9-1 Hinged-Hinged Beam on Elastic Supports 69 10-1 Analyzer Circuit for Transverse Beam Displacement at the nth Station, Including Transverse Shear 75 10-2 Comparison of Transverse Shear Effect on NormalMode Frequency for a Uniform Free-Free Beam 76 A-1 Comparison of Mode Shapes for Cellular and Continuous Clamped-Clamped Beams A-17 12

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX II Outline for Technical Report AIR-8 Entitled ELECTRONIC DIFFERENTIAL ANALYZER SOLUTION OF BEAMS WITH NONLINEAR DAMPING TABLE OF CONTENTS Chapter Title Page PREFACE ii ILLUSTRATIONS v 1. INTRODUCTION 1 1. 1 Equation for Lateral Vibration of Beams 1 1. 2 Finite Difference Method for Approximating Derivatives 5 1. 3 Princiiples. of Operation of the Electronic Differential Analyzer 5 2 CANTILEVER BEAM WITH.VELOCITY-SQUARED DAM1PING 10 2. 1 Beam Equation Including Velocity-Squared Damping 10 2. 2 Equivalence of Damping-Coefficient Size and Amplitude of Vibration 10 2. 3 Difference Equations for the Cantilever Beam with Velocity-Squared Damping 11 2. 4 Analyzer Circuit for the Cantilever Beam with Velocity-Squared Damping L3 2. 5 Damped First-Mode Oscillation 14 2. 6 Approximate Theoretical Solution 14 2. 7 Impulse Response of the Cantilever Beam with Velocity-Squared Damping 17 3 CANTILEVER BEAM WITH COULOMB DAMPING 3. 1 Beam Equation Including Coulomb Damping 22 3. 2 Difference Equations for the Cantilever Beam with Coulomb Damping 23 3. 3 Analyzer Circuit for the Cantilever Beam with Coulomb Damping 23 3. 4 Impulse Response of the Cantilever Beam with Coulomb Damping 24 BIBLIOGRAPHY 13

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN ILLUSTRATIONS Figure No. Title Page 1-1 Cantilever Beam 3 1-2 Cantilever Beam Divided into Stations 6 1-3 Operational Amplifier 7 1-4 Servo Multiplier 8 2-1 Analyzer Circuit at the nth Station for the Cantilever Beam with Velocity-Squared Damping 13 2-2 Damped First-Mode Oscillations of Uniform Cantilever Beam with Velocity-Squared Damping 15 2-3 Variation of Logarithmic Decrement 8 with Amplitude of Oscillation 18 2-4 Unit Impulse Response of 5-Cell Uniform Cantilever Beam with Velocity-Squared Damping; Displacements at Stations 2, 3,4, and 5 20 2-5 Unit Impulse Response of 5-Cell Uniform Cantilever Beam with Velocity-Squared Damping; BendingMoment at Stations 1,2, 3, and 4 21 3-1 Analyzer Circuit at the nth Station for Cantilever Beam with Coulomb Damping 24 3-2 Impulse Response of a 5-Cell Uniform Cantilever Beam with Coulomb Damping 25 14

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX III Outline for Technical Report AIR-10 Entitled APPLICATION OF DIFFERENCE TECHNIQUES TO HEAT FLOW PROBLEMS USING THE ELECTRONIC DIFFERENTIAL ANALYZER TABLE OF CONTENTS Chapter Title Page PREFACE 1 INTRODUCTION 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 3. 4 Solution in a Two-Dimensional Homogeneous Medium 38 15

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 4 HEAT FLOW IN CYLINDERS 43 4. 1 The Heat Equation in Cylindrical Coordinates 43 4. 2 Solution by Separation of Variables 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 Spherically-Symmetric 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 SEMIINFINITE MEDIA 80 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 DifferenceEquation Solution 90........ 16

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 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 17

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN ILLUSTRATIONS Figure No. Title Page 1-1 Temperature Distribution 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 /t for One-Dimensional Heat Flow n 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 Initial Temperature 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 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 18

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN 4-2 Comparison of 4 - 1/2 - Cell Mode Shapes with J, Modes 1 and 2 52 4-3 Comparison of 4 - 1/2-Cell Mode Shapes with J, 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 u./2(0) = 1 56 4-6 Axially -Symmetric Temperature Distribution with 11/2(0) = 1 57 4-7 Axially-Symmetric Temperature Distribution with U2-1/2(0)= 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)= 1 60 4-10 Axially-Symmetric Temperature Distribution with 5-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 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-Cymmetric 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 = x- 75 6-2 Analyzer Solution for Heat Flow in a Uniform Slab; Equal Stations along y, where y = fE 77 19

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN 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 U = 0.2 92 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 for Checking the Exact Solution by Separation of Variables 96.... — - 20

BIBLIOGRAPHY 1. D. W. Hagelbarger, C. E. Howe and R. M. Howe, "Investigation of the Utility of an Electronic Analog Computer in Engineering Problems," Ext. Memo. UMM-28, Univ. of Michigan Eng. Research Inst., AF Con. W33(038)ac-14222 (Project MX-794); April 1, 1949. 2. C. E. Howe, R. M. Howe and L. L. Rauch, "Application of the Electronic Differential Analyzer to the Oscillation of Beams, Including Shear and Rotary Inertia," Ext. Memo. UMM-67, Univ. of Michigan Eng. Research Inst; January 1951. 3. Howe, Howe and Rauch, "Solution of Linear Differential Equations with Time Varying Coefficients by the Electronic Differential Analyzer" Proj. Cyclone Symposium II, Part 2, 89; May, 1952. 4. G. M. Corcos, R. M. Howe, L. L. Rauch and J. R. Sellars, "Application of the Electronic Differential Analyzer to Eigenvalue Problems," Proj. Cyclone Symposium II, Part 2, 17; May, 1952. 5. R. M. Howe, "Application of Difference Techniques to Heat Flow Problems Using the Electronic Differential Analyzer" AIR-10, U.S. Army Contract No. DA-20-018-ORD-21811, Univ. of Michigan Eng. Research Inst; May, 1954. 6. C. E. Howe and R. M. Howe, "Application of Difference Techniques to the Lateral Vibration of Beams Using the Electronic Differential Analyzer," AIR-7, U. S. Army Contract No. DA-20-018-ORD-21811, University of Michigan Eng. Res. Inst; Feb. 1954. 7. R. M. Howe, "Electronic Differential Analyzer Solution of Beams with Nonlinear Damping," AIR-8, U. S. Army Contract No. DA-20-018-ORD21811, University of Michigan Eng. Res. Inst.; April, 1954.