SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS BY DIFF- ENCE:ETHODS USING THE ELTCTRONIC DIFFIRE`iNTIAL ALNALYZEI by V. o. HIanenman, Jr. Research Engineer,.Sngineering Research Institute University of Michigar t. 1., Howe Assistant Professor of Aeronautical -Engineering University' of li chigar Departlrient of Aeronautical Engineering Engineerini- g iesearch Institute University of lichigan AIR-t-1 October 1951

^AIwRl | ENGINEERING RESEARCH INSTITUTE A-1 _ _UNIVERSITY OF MICHIGAN P1 PREFACE The purpose of this report is to summarize the results of an investigation into the suitability of an electronic differential analyzer for solving linear differential difference equations. These equations arise when a partial differential equation is converted to a system of ordinary differential equations by replacing one or more of the partial derivatives by finite differences. The report includes both a theoretical analysis of the accuracies attainable using the difference method and actual examples of solutions of specific problems by the electronic differential analyzer. Three general type of partial differential equations are included; the heat equation, the wave equation, and the vibrating beam equation. No attempt is made to discuss in a detailed manner the theory of operation of the electronic differential analyzer, nor are the actual circuits of the d-c amplifiers and power supplies given here. For this information the 1,3 reader is referred to other reports. 3 The actual computer solutions were carried out on the electronic differential analyzer of the Department of Aeronautical Engineering. As a result of the promise shown by the difference method discussed in this report, construction of an 80-amplifier analyzer has begun. Complete details of this new computer will be presented in a later report.

A i | ENGINEERING RESEARCH INSTITUTE P a ______ AIR-1 v**>l _________ UNIVERSITY OF MICHIGAN Page TABLE OF CONTENTS Chapter Title Page 1 IINTRODUCTION 1 1.1 Usual Differential Analyzer Technique for Solving Partial Differential Equations 1 1.2 Replacement of Partial Derivatives by Finite Differences 1 2 PRINCIPLES OF OPERATION OF THE ELECTRONIC DIFFERENTIAL ANALYZER 4 2.1 Introduction to Operational Amplifiers 4 2.2 Solution of an Ordinary Differential Equation with Constant Coefficients 7 2.3 Solution of Two Simultaneous Ordinary Differential Equations 12 3 SOLUTION OF THE HEAT EQUATION 16 3.1 Equation to be Solved 16 3.2 Derivation of the Difference Equations 17 3.3 Imposing Boundary Conditions 18 3.4 Imposing Initial Conditions 20 3.5 Complete Differential Difference Equations for a Given Set of Boundary Conditions 20 3.6 Theoretical Solutions of the Difference Equations for Heat Flow 23 3.7 Computer Solution for One-Dimensional Heat Flow 32 3.8 Computer Solution for Two-Dimensional Heat Flow 36 3.9 Summary of Investigation of the Use of Difference Techniques for the Heat Equation 41 4 SOLUTION OF THE WAVE EQUATION 43 4.1 Equation to be Solved 43 4.2 The Wave Equation for the Stretched String 44 4.3 Derivation of the Difference Equation for the Stretched String 45 4.5 Solution of the String Equation by Separation of Variables 46 4.6 Difference Equation Solution of the Stretched String 49 4.7 Summary of Investigation of Difference Techniques Applied to the Wave Equation 50 5 VIBRATING BEAMS 52 5.1 Equation to be Solved 52 5.2 Derivation of the Difference Equation for the Vibrating Beam 53.'LILJII..II 1 1 I I I I I I. I I L.I.II

ENGINEERING RESEARCH INSTITUTE Pa ___ AIR-1 l"_ _UNIVERSITY OF MICHIGAN.g Chapter Title Page 5 (cont' d) 5.3 Representation of Boundary Conditions in the Difference Equations 54 5.4 Computer Circuit for Solving the Cantilever Beam by Difference Techniques 55 5.5 Theoretical Solution of the Difference Equation for Vibrating Beams 57 5.6 Computer Solution for a 5-Cell Cantilever Beam 67 5.7 Accuracy Requirements for the Difference Technique 72 5.8 Summary of Difference Technique for Vibrating Beams 74?,L

ENGINEERING RESEARCH INSTITUTE AIR-1 Page 1 UNIVERSITY OF MICHIGAN CHAPTER 1 INTRODUCTION 1.1 Usual Differential Analyzer Technique for Solving Partial Differential Equations The electronic differential analyzer is limited to the solution of ordinary differential equations. When one desires to solve a linear partial differential equation on the analyzer, it is necessary to separate variables and hence convert the partial differential equation to ordinary differential equations of the eigenvalue type. One must then find, usually by trial and error techniques, the normal modes, from which the complete solution to the original problem can be built up. The above method of separating variables and obtaining a series type of solution can be carried out fairly efficiently on an electronic differential analyzer. 1''3 Certainly, for most problems the analyzer is much faster than any hand methods. But for the engineer who is interested in getting numerical answers to specific problems, even the analyzer approach might seem somewhat tedious. It therefore would be highly advantageous to be able to solve the partial differential equations directly. This can be done by replacing some of the partial derivatives by finite differences in order to convert the original partial differential equation into a system of ordinary differential equations. In the next section we shall show how this is done. 1.2 Replacement of Partial Derivatives by Finite Differences 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 y1 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 -j! (i.e., the partial derivative of y with respect to x at the 1 station) 2is given by the difference i i i i i i..

ENGINEERING RESEARCH INSTITUTE AIBR-i |UNIVERSITY OF MICHIGAN 2 _ZI. ^- o. (1-1) ax 1/2 AX * In fact the limit of equation (1-1) as Ax -0 is just the definition of the partial derivative at that point. Writing (1-1) in more general terms ~s I,y-nynx n-1/2 x(1-2) In the same way'n2i/2 - I- ^~ ^3L | (.-3) ax n n+/2 x n-1/2 or from equation (2-2) a21 1- 12y +y - a --- 1. (1-4) 3x ( x)2 Similarly ( n Similarly _,IY_ Yn+l 3Yn + 3n- + n-2 I (Ax)3(1-5) 3 n-1/2 (x) and Ia Yn+2 -'4n+l+ 6Yn - Yn-l Yn-2 (ax4 (Ax)4

A|p Il ~ENGINEERING RESEARCH INSTITUTE l ___AIRi-1 ____ |__________ UNIVERSITY OF MICHIGAN Pag 3 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), yl(t),...yn(t),.... Before considering how a specific partial differential equation is transformed into differential difference equations or how boundary conditions are imposed, we will review the principles of operation of the electronic differential analyzer. The reader familiar with such a computer may choose to omit the following chapter and go immediately to Chapter 3 (Solution of the Heat Equation).

ENGINEERING RESEARCH INSTITUTE Pa UNIVERSITY OF MICHIGAN CHAPTER 2 PhIINCIPL~E OF OPihAT1ION uF Tii.J LiCTiuLiIC LIFi'ihtATIAL ANALYZIE 2.1 Introduction to perational _ plnifiers The basic component of the electronic differential analyzer is the operational amplifier, which is shoawn scheiatically in Figure 2-1. It consists of a d-c voltage amplifier of high gain, an input impedance Zi, and a feedback imtpedance Zf. i2 + zi — em ~ ~~e' D.C. AMPLIFIER e2 Figure 2-1. Operational Amplifier. If we neglect the current into the d-c amplifier itself (i.e., neglect the current to the grid of the input tube), it follows that i, = i2. Let us also neglect the voltage input e' to the d-c amplifier in comparison with the output voltage e2 or the input voltage el to the operational amplifie We then have iI = i2

ENGINEERING RESEARCH INSTITUTE P AIR-1 Page l _ A - IR-i_ _UNIVERSITY OF MICHIGAN or el e2 Zi Zf from which f e = - ~ e (2-1) which is the fundamental equation governing the behavior of the operational amplifier. In general, Zf/Li is made the order of magnitude of unity. We shall now consider the scheme by which the operational amplifier can be used to perform three different functions. (a) Multiplication by a constant. If we wish to multiply a certain voltage el by a constant factor k, we need only make Zf/Zi = k. From equation (2-1), then, the output voltage e2 of the operational amplifier will be given by e = k e. (2-2) Thus the required multiplication by a constant has been achieved, except for a reversal of sign. For example, if we wish k to be 10, we may let Z. = 1 megohm resistance, Zf = 10 megohms resistance. If we also desire the sign of e2 to be the same as el, we must feed e2 through an additional operational amplifier with 2. = Zf = 1 megohm. This second operational amplifier merely acts as a sign changer by multiplying any voltage by -1. (b) Addition. In order to add a number of voltages, say e, eb, and ec, the arrangement shown in Figure 2-2 is used. Here i + i + i = i2 and if we neglect e' as small compared vith e2, we have ea eb e e2 a Zb c f

ENGINEERING RESEARCH INSTITUTE Page ~~~AIR-1 ___ IUNIVERSITY OF MICHIGAN or zf Zf f e2 = [ +'b z e (2-3) a b c Thus the output voltage e2 is the sum of the three input voltages, each multiplied respectively by a constant -Zf/Zn (n = a, b, or c). The operational amplifier can, of course, be used in general to sumL any number of input voltages. + c —-- ZaZ Zf el D.C. AMPUFIER e2 +- _ Zc'1'l,, Figure 2-2. Operational Amplifier Used for Surixiation. (c) Integration. If we make the input impedance Z. a resistor and the feedback impedance Zf a capacitor, then the operational amplifier serves as an integrator. Referring to Figure 2-3, we see that if we neglect e' and let iI = i2 as before, we have e2 = -/ —- and i =-,. l lR

AIR-.1 ENGINEERING RESEARCH INSTITUTE Pa UNIVERSITY OF MICHIGAN from which e = - edt (2-4) 2= - 1C/ *< The output voltage e2 is then the integral with respect to time of the input 22 voltage e1 (multiplied by a constant factor -1/RC). R i +V e, e' D.C. AMPLIFIER e2 Figure 2-3. Operational Amplifier as an Integrater. 2.2 Solution of an Ordinary Differential Equation with Constant Coefficients In order to demonstrate how operational amplifiers performing the above three functions can be combined to solve ordinary linear differential equations, we will now set up the amplifier circuits required to solve the following differential equation: a2 *- al ldt aoy f(t) (2-5) a2dt 2 ao = f tdt subject to the initial conditions y(O) =V1 and A (0) =-V2. (2-6)

AIR-1 ENGINEERING RESEARCH INSTITUTE Page UNIVERSITY OF MICHIGAN The constants a2, al, and a are assumed positive. Since the electronic differential analyzer integrates with respect to tine, the independent variable t in equation (2-5) will be time. The dependent variable y is represented by voltage. The computer circuit for solving equation (2-5) subject to initial conditions (2-6) is shown schematically in Figure 2-4. If we assume that the output of amplifier A2 is a2y, then this voltage is multiplied by -1/a2 and integrated once in passing through amplifier A3, the output of which is therefore -y. This voltage is in turn multiplied by -1 and integrated once to give y as the output of A. In order to obtain +aly instead of -y, it is necessary to pass -y through summing amplifier Al. At the same time f(t) is fed into Al, so that the output of A1 is + aly - f(t). This output is then added to a y in amplifier A2, which finally has as its output - aly - aoy + f(t). But we originally assumed the output of A2 to be a2y. Hence a2y = - aly - aoy + f(t), which is just the equation which we wish to solve. 02 ~ rllI $~~-f(t) A-' aA W< I/ o f(t) ALL RESISTOR UNITS ARE MEGOHMS ALL CAPACITORS ARE IMFD. GROUND CONNECTIONS ARE OMITTED FOR CLARITY Figure 2-4. Computer Circuit for Solving a2y + aly + aoY = f(t).

AENGINEERING RESEARCH INSTITUTE Page l AIRH~-i _ l _ _ UNIVERSITY OF MICHIGAN The initial conditions (2-6) are imposed as voltages across the integrating condensers of A3 and A. in Figure 2-4. When the two switches holding the initial voltages across the condensers are opened simultaneously, the solution of the problem begins, i.e., the time variation in voltage output of A4 represents y(t). The mechanical analogy to equation (2-5) is the mass-spring-damper arrangement shown in Figure 2-5. The type of transient response of such a system depends upon the damping ratio i, which is defined as al i2 Ti;oa2 (2-7) For <C1 the transient response is oscillatory (underdamped), and for ~ > 1 the transient response is exponential (overdamped). In order to have a specific problem for the computer, let us examine the response of our second-order system to a step-function input force. For simplicity, assume a = a2 = 1. Then, = a1/2 and we can vary f merely by changing a1. Also, let us assume that initially y(O) = y(0) = 0. y~~~~ O//// ///,' 0 SPRING AMPING CONSTANT CONSTANT MASS OISPLACEMENT 02 y(t) APPLIED FOIGE fit) Figure 2-5. One-Degree-of-Freedom Mass, Spring, Damper System. I 1 i i iii~i ii ii

ATlil ENGINEERING RESEARCH INSTITUTE AIR-i I ____ UNIVERSITY OF MICHIGANge 10 CO.1 PRINTED IN USA S LV div/SeC __~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~_ ii 7* —I 0.2, 1 di'/sec ^g'-M ^-n^^^ -rf^ "*~~^~~~~~~~~~~~~~~~~~~r —-' — 4 ^:'^g'"~~j^ --'~~i~i~t' 14 -- Y^ ^-^*^ ^.~~ ~~ ~~-i'-!-t —- Ir ~ \ -I, —- ~ —--- i — ^^g^^^7?^^^^^/i-i__i / i /-r-i-ff-?''' —1H 4 1i-l\ iFigue 2 S R C0.1 PRINTED__..L! IN U i div ise ~~_l~~~ri _ ~ ~ _________111,!-1 —-T — T f ~ ~ ~ --- -` CHART NO. BL-924 THE BRUSH DEVELOPME -- i~~~~~~T __t ___~~~~~~~~~ \_: L\I:_ Figure 2 6. Step-Function Response ofl Second-Order System.

~AIR-i~ ENGINEERING RESEARCH INSTITUTE Page ~~~~AIR-1 U__ UNIVERSITY OF MICHIGAN The computer response to a step-voltage input at t = 0 for 3 = 0.2, 3 = 1, and = = 2 is shown in Figure 2-6. The input and output voltages y(t) were recorded on a Brush, Model BL-202, two-channel magnetic oscillograph. In Figure 2-7the response to a step-function input is shown when 4 = 0. The result is undamped simple harmonic motion of period 21T seconds. One can obtain some measure of the inherent accuracy of the electronic differential analyzer by observing that it takes many hundreds of cycles before this sinesoidal oscillation has decreased to half-amplitude. Fg 0 ur 21 Seo d- rv/ecryse wCHART NO. E.....1 m I 1! I Figure I2\ -. Seond-Orer -. i t =-\ \ i\ —-- \ D pi —.

ENGINEERING RESEARCH INSTITUTE _ AI1 ___R-i) |._ UNIVERSITY OF MICHIGAN 2.3 Solution of Two Simultaneous Ordinar Differential Equations Consider the two-degree-of-freedom system shown in Figure 2-8. The two masses ml and m2 are supported by springs having constants k and k2, and are coupled by a spring with constant k3. If the masses are constrained to vertical motion, then in terms of the deflections y1 and Y2 we find the following equations of motion: ml + k, + k )Y1- k3Y2 0 and (2-8) m2 + (k2 + k3)y2 - k3y = 0. The computer circuit for solving these equations is shown in Figure 2-9. The analysis of the circuit is exactly the same as the analysis described for Figure 2-4. The upper bank of amplifiers in Figure 2-8 represents the first equation of (2-8), the lower bank the second of (2-8). Cross-connections between the two banks are made as necessary to satisfy equation (2-8) For simplicity, assume ml = m2 = 1, k = k2 = 1. Let us make k somewhat less, say 0.2. For initial conditions we will set Yl(0) = y = 0 (i.e., no initial velocity for the masses). As a first example we let the masses start with equal displacement, so that Yl(0) = Y2(0) = Vo The response of the masses for these initial conditions is shown in Figure 2-10o Evidently the masses perform simple harmonic oscillation in phase with each other, so that the coupling spring is never compressed or expanded. The period of oscillation of each of the masses is therefore the natural period of the single mass-spring system without coupling, namely 21 seconds. This is evident in Figure 2-10. For a second example, let the masses start with equal but opposite displacement, so that Yl(0) - Y2(0) = V. The response for these initial conditions is shown in Figure 2-11, and now the masses oscillate 180 degrees out of phase and with a shorter period (actually 2TT/g/24 seconds). In the above two examples we have found the two normal modes of vibration of our two-degree-of-freedom system. Any other motion which the two masses can exhibit must consist of a superposition of these two normal modes. For example, if we start one mass with finite displacement and the other mass with zero displacement, so that Yl(0) = V, Y2(0) = O, we get the response l

A.-1 -ENGINEERING RESEARCH INSTITUTE Page 13 AlHR-1 ___ UNIVERSITY OF MICHIGANPae 3. k! m yj y? k3 m"2 - 2 k 2 Figure 2-8. Two-Degree-of-Freedom 3ystem with 3pring Coupling. k, + k3 V I~~~~~~~~~-V k3 ml m -h ~ ~ ~ ~~ I -y yj k2+ k3: lw l w w k f -^~~~~~~~~s MAIRA-~ Figure 2-9. Computer for 3olvinc the.jir;raltaneous Differential equations of 1Bqua.tion (2-8).

ENGINEERING RESEARCH INSTITUTE p AIR-1 | UNIVERSITY OF MICHIGAN Page 14 J I I I I I -i f- / 1 1 1 1W 1 / I t t I I I XII 1 1 1 1 1 1 1 0PRINTED.IN U.S.A. d sec Figure 2-10. Coupled System, Iasses Started with Same Deflection THE BRUSH DEVELOPMENT CO. PRINTED IN U.S.A. / i X r X X - t- I.i Figure~~~~~~~~~~~~ 2-1.ope ytm lse tafe ihOpst elcin

AIR-1 ENGINEERING RESEARCH INSTITUTE Page _ __ AIR-i ___ _UNIVERSITY OF MICHIGAN 15 shown in Figure 2-12. Evidently, the masses are exchanging energy back and forth through the coupling spring. The frequency of energy exchange is just the beat frequency of the two normal modes. It is obvious that an arbitrary forcing function could be added to the two masses in Figure 2-8 along with damping. Such effects can readily be incorporated into the computer circuit exactly as done in Figure 2-4. This simple example of a two-degree-of-freedom system illustrates the way in which the differential analyzer can be used to solve simultaneous ordinary differential equatisns and hence differential difference equations. We now turn to a consideration of a sirnple partial differential equation, namely, the heat equation, which will be solved by the difference technique. CHART NO. BL 909 THE BRUSH DEVELOPMENT CO. PRINTED I Figure 2-12. Coupled System, One Mass Started with Zero Deflection. 7\\B A A \ \ \ \ \ \ \ \ \ \t \ \ \ \ \ \ \ \ \ \ \ \ \CHART~~~~~ NO L99 TEBUHDVEOMN O RNE I I / I I I J I I I / I I L ff I I I I I I

ENGINEERING RESEARCH INSTITUTE 1 P1' ________ UNIVERSITY OF MICHIGAN Page 16 CHIlFTIFR 5 SOLUTION (OF THE THAT EQUATION 5.1 Equ3tion to be Solved As a first example of a partial differential equation, it seemed advisable to select the equttion of heat flow through a continuous medium, since it involves second order spacial derivatives and only first order time derivatives. The basic heat equation is given by cl ( V- KrVu + f (3-1) 2) where u = temperature and is a function of the spacial coordinates and time T = thermal conductivity, in genera1l a function of the spac ial coordinates, C = specific heat, a function of spacial coordinates, = density, also a function of spacial coordinates, t = time, f = rate of heat supplied by sources in the medium, a function of spacial coordinates and time. The left-hand side of equation (5-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 7. K7u term) and second due to the heat flow into the volume from sources within the volume itself (the f term). Written in terms of Cartesian coordinates x, y and z, equation (3-1) becomes KKK (5-2) dt ax ax ay ay az az The a ctual heat flow or flux due to conduction normal to any unit surface is given by - K X (component of V u normal to the surface). Thus the heat flux F^ across a unit surface normal to the x direction is given by

l A lR- ENGINEERING RESEARCH INSTITUTE P 17 _AIhRn-1 __1 _UNIVERSITY OF MICHIGAN Page F K _; (3-3) X ax In a given heat flow problem it is necessary to stipulate spacial boundary conditions either on the temperature u or the heat flow -K7u, as well as the initial temperature distribution throughout the medium. 3.2 Derivation of the Difference Equations For simplicity in illustrating the application of difference techniques, let ts assume for the time being that spacial variations in the temperature u are confined to the x direction. Equation (3-2) then becomes C((x2)1u(xt) -a u(xt) + f(x,t), (3-4) clt =x ( x where we have let C(x) = c(x) (x). (3-5) For example, this could represent the temperature distribution in a medium between two infinite slabs, as shown in Figure 3-1. Following the technique discussed in Section 1.2, we will consider only values of u at certain equally spaced stations along the x coordinate axis. Thus u(x,t) is replaced by ul(t), u2(t),...., etc. If Ax is the distance between stations, we can write for the heat flux Fnj at the n-~ station CONTINUOUS MEDIUM CONDUCTIVITY ~ K() / HEAT CAPACITY* C(l) TEMPERATURE * att) Figure 3-1. Temperature Distribution Between Two Infinite Slabs.

ENGINEERING RESEARCH INSTITUTE Page 18 _____ AIRtM"-1 ____________ UNIVERSITY OF MICHIGAN _ F 1 -K- - I -(U- n U l) n (3-6) n- ax x n n-l (un-3-6) In the same way, (K - ) at the nth station becomes a9 x <x axaax K n+ | K n | (Kau ax x a x ^ Ax or K n+ ( K na ( au, U n+l Un)- ax n Un-). (3 7) 0x axA From equations (3-4) and (3-7) we can now write the equation of heat-flow balance at the nth station. Thus du K1 K c dn 2,, (u- + (3-8) n dt (A X)2 n+l n ( x)2 n n(3 n where Cn is the heat capacity at the nth station and fn is the rate of heat supplied by a heat source at the nth stationl(f will in general be a function of time). Note that du /dt is now a total derivative and not a partial derivative, since by definition x remains fixed while we take dun/dt. Equation (3-8) will be iterated for different values of n until the boundaries in x are reached, at which point it is necessary to impose boundary conditions. 3.3 Imposing Boundary Conditions (a) Conditions on the temperature. Suppose that one of the boundary conditions specifies the temperature at x =0 (i.e., at the zero station). Then we have us = constant and hence F -x [u(t) -u (3-9)i

~AIR~ 1 ENGINEERING RESEARCH INSTITUTE P.____AIR-i7 __ E UNIVERSITY OF MICHIGAN P 1 Our difference equation for the first station is now du1(t) _ K Fu(t) u(t)) 1C 32 iS (U ul(t)?' T^ - l-(t)- u + fl(t) (3-10) C1 dt (a )2 2(t) - (t x)2 1 The equation for the second station proceeds as in equation (3-8) for n = 2. All we have done in imposing the boundary condition, then, is to fix u (t) at a constant value of u. If the temperature is specified at x = L (i.e., at the Nth station, where N = L/Ax), then for uN(t) we substitute uN = constant, the desired temperature. (b) Conditions on the heat flow. Often a condition is placed on the rate of heat flowing past a boundary, either that this flow be zero (as for an insulating boundary) or a constant. Supose we let F., = constant. Then the equation for the first station is du1(t) F3/2(t) - Fl C - 2 f(t) (3-11) 1 dt x +1 or dul(t) K F u(t)1 Fi C1 dt (x2 2(t) - ul(t) + + fl(t)' (3-12) The equations for u2, u3,... are the same as usual. If we desire FN+ = constant as a boundary condition, then the equation for the Nth station becomes duN(t) FN —~ K ut N dt A x t2 f (t)( u-lj The process of setting in boundary conditions is evidently quite straightforward. Notice, however, that when we denote temperature at integral I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ENGINEERING RESEARCH INSTITUTEp _____ AIR-1 ___UNIVERSITY OF MICHIGAN a stations, the boundary occurs at an integral station when temperature at the boundary is specified, whereas the boundary occurs at a half-integral station when the heat flow at the boundary is specified. 3.4 Imposing Initial Conditions In addition to specifying boundary conditions in tnis type of heat problem, it is necessary to specify the initial temperature distribution in our medium. Thus we have u1(0) = U1 u(0) = U2 u^(0) = U 3 3 (3-14) uN(O) = UN These initial conditions must then be imposed on the electronic differential analyzer, similar to the way in which initial conditions are applied in Figure 2-4. 3.5 Complete Differential Difference Equations for a Given Set of Boundar Qonditions For purposes of illustration, let us assume that the boundary conditions of our conducting slab in Figure 3-1 are that at x = 0 the tempera ture-remains fixed at u, and at x = L = A x(N+i-) the heat flow is zero. The space between x = 0 and x = L is therefore broken into N cells, and from equations (3-8), (3-10), and (3-13) we have for our complete set of differential difference equations du K K c d = (^x)2 2 (u -2 f 2 dt ( x)2 (U3 -u) (x2 (U2 ) + f du2* _ -

l A-IR 1,| ENGINEERING RESEARCH INSTITUTE Page 2 _____<u'"L ______AIR-i 1 | ___ [UNIVERSITY OF MICHIGAN * * * 0 du K 1 K 1 C - n —xa- (u) - 2 ( -u Cn dt 2 (A ) 2 n+l n (Ax)2 n n-) n ~~* * ~~~~~.... (3-15) * * * * * * 0 -KN - K /) -3/ 2 - f) - -1 d-dt (Ax)2 -1 (Ax)2 n-1 n-2 N-2 duN K 1 - 2 —(u -u )+f CN dt (x)2 (Un n-) + The initial conditions specify the temperature for each station at t = 0 &Iee equation (3-14)]. A schematic diagram showing all the locations relative to the conducting slab is shown in Figure 3-2 for N = 10. / / \SU^S^I^^^^u^^: INSULATOR < T O z L ll Figure 3-2. Station Arrangement for N =10. In Figure 3-3 the computer arrangement for solving equations (3-15) is shown. Note that the outputs of each successive row of amplifiers are Rao~~~ ROL

AI,-p~i,ENGINEERING RESEARCH INSTITUTE __ AnIR ________ UNIVERSITY OF MICHIGAN g _ K,,2 Ax ____^ 1^t/AX Ax L~Ul^~ ~~ K3/ Ax i\~~~~ -i fg -C4u, V~~\ y r Ax * 1^ ~Ax ^ - r'/i+ FMFiu.Com C1 ve eel 1a with Teperatur K 5/2- adetFlx=0Atx + L+/)x A x JX/d C3 -In-~~~~~~~Fn / - rr -4P U'\ r^^ " /( i I II I I~~~~A f v/ith T = 0 1t 0 e0( CN-1 6N-1 - -FN x + fu_3/ FNxi GM-Fn /

ENGINEERING RESEARCH INSTITUTE Page.__A-ir-l _____ _UNIVERSITY OF MICHIGAN IP reversed. This allows the necessary differences to be tak-n:J thout signreversing amplifiers. Note also that the heat flow or flux F is available at any half-station as a dopend.laent variable. Thus one can observe directly as a function of t tihe te temperature u and heat flux F across the slab. It is po3sible to reduce the number of al;plifiers shown in Figure 3-3 from three to one per station. The exact way in which this is done is explained in Section 3.7. In many ways, however, the circuit of Figure 3-3 is simpler dx:spite the increased nurib:r of amplifiers. To change the conductivity K or heat capacity C at any station, one need only vary the appropriate resistor. Initial terperature distribution across the slab is changed by settin the U1, U2,... U voltages to the desired values. The heat sources through the slab are represented by the voltages fl, f2,.. fI,i'hiih Tray be varied a- a function of tie in any aesired manner. Before actuCilly solving a p -.'rticular heat problem with the differential analyze-, it nmiht be.rell for uis to;aie soxmle calculations of the accuracies we can oexpect wi;n using thze &iff.- ence technique. 3.6 Theoretical S3lutions of the D iffrence E curations for HIeat Flow (a) Preliridnary solution by separation of variables. In order to evaluate the accuracy of the difference technique, it is -worthwhile lfirst to solve the partial differential eju-lations of heat flow by sep-ratinj variableso For - i;.plicity we will solve the probleirn of the teimperature distribution bet;ee~ -ltzo irniite s.laLs hllc a't a terr.erat.re of zero (Figure 3-4). As.ume that th~ me diuui has constant corductivity;-:,nr constant sipccific h:^at capacity c. Also assule no heat sources within the medium. Then from equation (3-4) cs au a2 ( - i ^' 2 The boundary conditions are u(O,t) = u(L,t) = 0. (3-17) Let us aseuime as a siipl e initiail condit ion that; the te.:iperature in the r:;ediuml is everywhere constant at t = O. Thus u(Qi,0) = U'-c lnst'i;:it.Lo (3-18)

ENGINEERING RESEARCH INSTITUTE Pa AIR-1i ___ I ___ UNIVERSITY OF MICHIGAN ag 24 KC- CONSTANT C: CONSTANT Tr. / / T x= O x L Figure 3-4. Temperature Problem Solved by Separation of Variables. To solve this problem by separation of variables we let u(x,t) = X(x)T(t). (3-19) Substituting (3-19) in equation (3-16) we find that xlitW cs T__t) (3-20) X(x) K T(t) X KT X(xt) + o2X(x) =-0 For equation (3-23) to hold for all values of x and t, it necessarily follows 2 that both sides of the equation equal a constant, say -o. Thus X c eT — o (3-21) from which we obtain the following two equations: d (x) + cA(x) = o and (3-22) T'(t) K. o<(t) =0. I -_ I~~~~~c

- ENGINEERING RESEARCH INSTITUTE P EAN-1GUNIVERSITY OF MICHIGAN g The general solution to equation (3-24) is X(x) = A cos cKx + B sin o(x (3-23) while the solution to equation (3-25) is K 2t T(t) = De c. (3-24) The boundary conditions (3-17) can be satisfied only if X(O) = X(L) = 0. Reference to equation (3-23) shows that this is true when A = 0 and when ( has discrete values o( given by rIT n -L n =1, 2, 3. (3-25 These values o( are known as eigenvalues, and the resulting functions n Xn(x) = Bn sin nLx (3-26) are known as eigen functions or normal modes. They are orthogonal, since L I X(x)X( (x)dx = 0 n m (3-27) n- n=m. L = 2 n =m. From equation (3-19), (3-24), and (3-26) we see that the complete solution to our heat problems can be written as ~00 K (IK n)2 nu(xt) E= Bsin L e". (3-28) n=l

ENGINEERING RESEARCH INSTITUTE Pa 26 l AIR-1i EG NUNIVERSITY OF MICHIGAN To evaluate the B constants we must apply the initial condition (3-18). Thus 00oo niTx u(O,t) = U = Bn sin L. (3-29) n=l Multiplying both sides of equation (3-29 by sin m x and integrating between 0 and L we have L L Usina1 dx B sin - sin a -dx LL n L f U sin L xf s sin d. 0 n=l 0 From the orthogonality relation (3-27) it is evident that the right side of the above equation vanishes except when m = n, so that L U sin x dx = B 0 from which B = MU n odd (3-30) n nil = O n even The final solution can now be written in the series form. u(x,t) A 2 1 sin (2^j-1) x e c L t (3-31) j =1 This solution actually represents an infinite number of sinesoidal temperature distributions across the medium from x = 0 to x = L. At t = 0 the sine waves all add up to give the initial flat temperature distribution. For t >0 i - -. ilI_ II I I I 1 I I..,E. l ~ I l ilI I ~ I II I II I I -.. l

ENGINEERING RESEARCH INSTITUTE P UNIVERSITY OF MICHIGAN the sine waves decay exponentially at different rates, with the decay rate faster for those sine waves having more nodes and loops. The resulting temperature distribution at various times is plotted in Figure 5-5, where a dimensionless time variable 7has been used. 7' is defined as 7"' K t (3-32) cSL 2 Thus Figure 5-5 is independent of the physical constants of the problem. We will now proceed to calculate thed- Is for a differential n difference equation representation of the heat equations. If these( i's n agree well with the valuesO(n = nT/ L from the solution above, and if the equivalent normal modes Xn show good agreement with sine waves, then we can expect accurate results using the difference technique. (b) Solution of the difference equation for N cells. When the space between x = 0 and x = L in Figure 5-3 is broken up into N cells so that there are N - + 1 temperature stations, the general difference equation is given by (5-8). At station 1 and station N-l the difference equation is obtained from equation (5-8) by setting u and uN equal to zero respectively. In the problem under consideration the conductivity K and specific heat capacity C are constant. By proper choice of our distance variable x we can make.4 x = 1, so that L = NAx = N. By proper choice of our time variable t we can make C/K = 1 so that for f = 0 equation (5-8) becomes for the ith cell. u. = u+l 2ui -1 ui-l (3-3355) For N cells and for boundary conditions u = uN= -0, we hayve = Up - 2uU2 = u5 - 2u2 + u * **

T =0 O 1.0......... H.. H." I ^ ^ — ~"]fc^ ---- -.~~^ —^.^ ____ ~~~CONTINUOUS m 1.0 l |^~~~~~ T:~ ~ =0.01.'^ ^o I|^ ^MEDIUM, lT |'02/ t O. O < > H04 ^^ss^ *^ 0 SOLUTION, ___ is ______ /\ / ( ^_____________ ^S_\. \ \ \______ I-I x1134578 "LO.0 o; C+ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~, 1, 11,1,

1ENGINEERING RESEARCH INSTITUTE P ____AIR-1i__EUNIVERSITY OF MICHIGAN Pae 2 (3-34) N-2= UN - 2N-2 + N-3 N-1 = - 2uNl + %-2 To solve for the normal modes we assume that the ith temperature u. -At varies with time as a.e, where a. is a constant. If this is true, then equations (3-34) become (-+ 2)a1 - a2 = 0 -a + (-A + 2)a2 - a3 = 0 -a- + (-A + 2)aN2 - aN-3 = 0 (-X 2)aaN - N= 0 The only nontrivial solution of equations (3-35) is obtained when the determinant of the coefficient vanishes. Thus (\-2) 1 0 1 ( -2) 1 0 * 1 (A-2) 1 0 * 0. 0 0 0 0 0 * * 0 0 0 0 0 0 = 0 * 0 0 0 0 0 0 0 1 (A-2) 1 1 ()^-2) 1 * * * 0 0 * 1 (%-2)

ENGINEERING RESEARCH INSTITUTE Pge lAIR-1 _______ UNIVERSITY OF MICHIGAN 30 This determinant, when expanded, becomes a polynomial in of order N-1. The polynomial will have N-l positive roots X which are the eigenvalues for our 2 N cell system, corresponding to the eigenvalues Xon of the continuous system [see equation (3-25)]. To solve this determinant for a specific N is very tedious, and to solve it in general would be next to impossible. The roots Xn can be found much more easily by the following procedure: Assume that the spacial mode shape for the difference equations is the same as for the continuous equation, ioeo, sinesoidalo If this is true, then for the temperature u. at the ith station we have n TTi n t u. = a sin e.(3-36) From simple trignometry it follows that nTli niT ui+ + u = 2a sin - cos n. (3-37) From equations (3-33) we have for the ith station nri n n niTi nTT n (2 A) sin - e = 2 sin -- cos -- e from which = 2(1 - cos ) (3-38) By substituting equation (3-36) into the first and last of equations (3-34) it is easy to show that the boundary conditions are satisfied. Thus, our assumed solution (3-36) is the exact solution, where the eigenvalues Xn are given by (3-38) Expanding equation (3-38) in a power series, we have n 21 (325) for since here L = N 12 = N (We assumed earlier thatx = 1) In the limit of infinitely many cells N equation (3-39) reduces to equation (3-25) for can, since here L = N Ax = N (lie assumed earlier thatax = 1)o i l lIII I

AIR-1 ENGINEERING RESEARCH INSTITUTE Pae 31 UNIVERSITY OF MICHIGAN The equivalent of the decay constant O2 is n In Figure 3-6 the 2 fl ~~n sn nn percentage deviation in OCn due to the difference method as a function of the number of cells N is shown. Note that the lower modes (lower values of n) require fewer cells to give accurate decay constants. NUMBER OF CELLS 2 4 6 8 10 12 14 16 18 20 22 24'b a a i 0 th f zr 4 6 I z 10 IL Z I0 0iue?6'ret Iv~to nteLcg Crs tQn 4 ~ ~ ~~~,~ Zl'nto fteNr~e fCls

A.yph-l I ENGINEERING RESEARCH INSTITUTE Page 2 I ____ UNIVERSITY OF MICHIGAN To suwllarize, we see that iiien the spacial d ivtieves, of te heat ecqu&tion are replaced by.iitC difference, the resultin noriai' r l- O shapes agree exactly, whereas tihe decay constants (eigenvalucs) for each node are so-ewliat smaller. This xcans that the higher modes will decay soinehat slowe whicn the differential difference ecuation approximiation is; dsedC. The error is bigger for hiither modes, but fortunately the higher uods are ogenerall —y muchl less important. 3.7 Computer Solution for One-Dimensional Heat Flow We now proceed to the computer solution of the one-L-dinm;ensional Iheat flowv problem considered in tie last saection, namely the teAperature distribution between two infinite slabs a distance L apart and.ith boudaries held - at zero teniecrature (See Figure 3-4). -i c' can. select the uistancc variable so that 6 x = 1 and hence L = ilAx = i h, hore &i is the numoer of cells. xAfter proper choice of the units of tiime t so tihat K/c S = 1, the basic heat equation becomes from (3-16) a 2 — ^ (3-40) at Zx which in terms of a difference equation is du n dt Unl -2 + U1i-l-41) For the problemi:n 3`cti._n 13.b the initial t;ipArature cistmIbution was a constant iU. Thus, Ie inve the initial conilditions u (0) = U = constant. (3-42) Boundary conditions are u = u = (3-43) Let us solve this heat problem with the differential analyser for 9 cells, From equations (3-41) and (3-43) the difference equations become

ENGINEERING RESEARCH INSTITUTE _AIR-1 ENGINEUNIVERSITY OF MICHIGAN Pag u1 = u2 - 2u1 u3 = U - 2u3 + u2 (3-44) * * * * U6 =7 - 2u6 + U5 u7 = U - 2u + u6 u8 =-2ug + u7 where dun/dt has been abbreviated as u. Equation (3-44) is subject to the initial condition (3-42). The computer circuit used to solve equation (3-44) is shown in Figure (3-7). Note that only one amplifier per cell is needed. The results of the computer solutions are shown in Figure 3-8, where temperature at stations 1, 2, 3, and 4 appears as a function of time. Since our initial temperature distribution is symmetrical about the station 4-1/2, as are our boundary conditions, the temperature distribution remains symetrical as a function of time. Thus u4 = U u3 =2 = u7, and u1= u8. In order to compare the computer results with the solution shown in Figure 3-5 for a continuous medium, we must convert our computer time units to the dimensionless units of Figure 3-5. Remembering that we chose computer time units so that K/cs = 1, we have from equation (3-32) L2 or since L = NAx where Ax = 1, T = 1 t. (3-45) N2.

UNIVERSITY OF MICHIGAN P 0.5 INITIAL CONDITION CIRCUITS OMITTED FOR CLARITY 0.5 -U6 0.5 I 0.5 U7 0.5 05I -I m -Us Figure 3-7. Computer Circuit for Solving 9-Cell ileat Problem

AIR-1 ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN Page 3 _- -AB^tE^.^^ \-V^^;H DEVELOPMENT CO. PRINTED IN U.S.A.. t!!_'__ t - t _ /___ 1 _ 1.... -I........ _ ==:- - _ _ - /; —:::: --::;- - -: —;/ --:- = —- / —-—: - ---- — T —.._,.,: \^\^~~~~-. \ \ \I<A,, -. ___ -., _ -._.. = - i 1 1 77,. -_ _ _: _ CHART NO. BL-924 THE BRUSH DEVEI.......................:... — _-::': _ —-~ —_- _ f. - t -. -' -- _ —-?_~~~~~~~~-_....-._....:. —::-..-_-..........t — t — t-7 t —,,,, 1 Figure 3 -8 Computer Solutio n for 9-Cell Heat Problem.

ENGINEERING RESEARCH INSTITUTE 36 AIR-i I Page _____AIR-l"1 _UNIVERSITY OF MICHIGAN Thus for our 9-cell problem we divide computer time t in Figure 3-8 by 81 to obtain the dimensionless time T of Figure 3-5. In this way points from the computer solution are compared in Figure 3-5 with the theoretical solution for a continuous medium. The correlation is evidently quite good, as we could have predicted from our theoretical work in Section 3.6. It has already been pointed out that the temperature distribution is symmetrical for this problem. Therefore, the heat flow will be zero at station 4-1/2, and the appropriate boundary condition can be established thereo If this is done, it is only necessary to solve the problem half-way across the distance between the slabs, the solution for the other half being symmetricalo In this case our difference equations become e u1 =u2 -2u1 2 = u3 - 22 + u1 (3-46) u3 = - 23 + u2 u = -4 + 3 subject again to initial conditions (3-42). The four-amplifier circuit required to solve equation (3-46) is shown in Figure (3-9). In the same way, if the initial temperature distribution in our homogeneous medium had been antisymmetrical with respect to station 4-1/2, we could have treated the problem for N cells by setting uo = uN/2 = 0 and solving the N/2-cell problem. Here we must obviously have an even number of cells to begin with, whereas in the symmetrical case we needed an odd number of cells. It is evident that by considering symmetry effects the number of amplifiers needed may often be cut in half. Furthermore, any arbitrary initial temperature distribution can always be split into a symmetrical and antisymmetrical formo The solution for each of these initial distributions can then be found, and since the equations are linear, the final solution is the sum of the two solutions. Of course this procedure will only work when the conductivity K and the specific heat capacity c S for the medium are

ENGINEERING RESEARCH INSTITUTE ____~AIR-1 __ _UNIVERSITY OF MICHIGAN 0.5 UI 05INITIAL CONDITION CIRGUITS I OMITTED FOR CLARITY _U2 0.5 Figure 3-9. Comnputer Circuit for Solvin- )-Cell Sy-..ja-etrical Heat Problem, constant (as in our problem) or sy-iniaetrical about the center of the medium. Also, the boundary conditions must be symmetrical. 3.8 Comuter Solution for Two-Dim.ensional Heat Flow It was felt that it wvould be interesting to solve at least one second order space problem using the difference technique.,jeccause only one amplifier per cell is required, the heat equation wx-as selected as the simplest two-dimensional problem for the electronic differential analyzer. Consider the homogeneous mediurm'of rectangular cross section shiolnl in`igure 3-10. Let the temperature u be a function of x and y and independent of the height Z. The walls are held at zero temperature and separated by X and Y respectively..Je can select thle units of time so that h/cs = 1. Thius, the heat ecuation becomes 2 2 au 92u u 9;t p2 + 2 (3-47)

AIR-1 ENGINEERING RESEARCH INSTITUTE Page 8 UNIVERSITY OF MICHIGAN /////// —/-,-/// I /7 — ri-a-c9 -i-i — U21 212 u23 /ith I budr t-nitins"~< I I I I / I u(-X /2,y~t = u "2=1 =, ",t3 = O l I /, /, /,, —-+ -+ —t-t —— / / X L.Z_! 1 3 I'/ / /// / -/-////// Figure 3-10. Two-Dimensional Problem in Heat Flow. with boundary conditions u(-X/2,y,t) = u(X/2,y,t) = u(x,-Y/2,t) = u(x,Y/2,t) = 0. The division of the medium into cells and the numbering system for the cells is shown in Figure 3-10. By selecting the units in x and y properly we can make Ax = Ay = 1. Thus for an N X N cell division X = Y = N. Since the initial distribution is symmetrical with respect to the x and y axes, we can solve the problem in one quadrant (See Figure 3-10), subject| to the new boundary conditions (- x,y,t) = u(x,,t) = 0 (3-49) u. (0,y,t) = (x,O,t) =0. (3-50) Ax i~~~~~~~~~~~~~~~

AIR-1 ENGINEERING RESEARCH INSTITUTE Page UNIVERSITY OF MICHIGAN ge 0.25'1 0.25 0.333 I,, ull I'-U12:' zU13 0.25 e,0.25: 02. 22 U23 0.333 0.333 0.5333 0.: 1 I -u321 0U33 Figure 3-11. Computer Circuit for Two-3Dimensional Heat Problem.

AIR-i i I ENGINEERING RESEARCH INSTITUTE 40 ___ AIR-1i__ [ __ _ UNIVERSITY OF MICHIGAN The difference equation for the temperature u becomes n,m du n- S 3n U - 2u + +u -2u + dt n+l,m n,m n-l,m n,m+l n,m n,m-l or du n.m Uu -4u + u + u dt Un+l,m nl,m nm n,m+l n,m-l (3-51) Let us consider the 7-cell by 7-~cell system shown in Figure 3-10. The boundary conditions of equation (3-45) are imposed by setting Uo0 - o02 = uo3 -uO = 20 O 30= ~ and the boundary conditions of equation (3-50) are imposed by setting %3 = u14, u23 U24' u33 = 34 u31 = u4l, u32 = u42, U33 U43 The complete set of 9 difference equations representing the upper left-band quadrant of Figure 3-10 then reduces to = u21 - 4Ul + U12 u12 = u22 - 12 + u11 + u13 13 U= 23- 3U13 + 12 121 U31 - 4u21 + U11 + 22 (3-52) U22 U32 + U12- u. 22 + U23 + 21 123'U3 +U 3 - 3 23 + 22 31 = 21- 3%1 + "32 u32 * 22 - 3U32 + u3 1 31 33 "32 - 2 u3 + 123.

ENGINEERING RESEARCH INSTITUTE.AIR-1 Page 41 _ AI 1 |______ UNIVERSITY OF MICHIGAN As in the case of the one-dimensional heat flow treated in the previous sections, we solve the problem for an initial temperature which is everywhere constant. Thus, the initial condition becomes u = U = constant. n,m The computer circuit used to solve the set of difference equations (3-52) is shown in Figure 3-11. Computer solutions for temperature ull, ul2, Ul3, and u33 as a function of time are shown in Figure 3-12. By symmetry U12 = u21, and ul3 = u31o The resulting temperature agrees closely with the theoretical distribution calculated for a continuous medium. 3.9 Summary of Investigation of the Use of Difference Techniques for the Heat Equation We have shown that it is both simple and straightforward to solve the heat equation with the electronic differential analyzer by replacing spacial derivatives with finite differences. Normal mode shapes show exact agreement with those calculated by separation of variables for the simple problems considered. Decay constants corresponding to the various modes also show good agreement but tend to be somewhat lower than the values calculated by separation of variables, particularly for higher modes or if fewrer cells are used. For most engineering problems the order of eight to sixteen cells per spacial dimension should be completely adequate (see Figure 3-5). Only one operational amplifier is needed per cell, although in some problems it may be more convenient to use three amplifiers per cell. The problem is completely stable, and the final outputs of the computer are temperature and heat flow as a function of spacial coordinates and time,

ENGINEERING RESEARCH INSTITUTE __ AIR-)1 I___ __ UNIVERSITY OF MICHIGAN age 2' _ _-'-'- 1 -- --- t i d _-I iART NO. BL-924 THE BRUSH DEVELOPMENT i-i I!< T i -j I I!, I F4i / I /': i'^ " j,; i' fI \\ \ ^ -+ — -------- ----— i. _ i. I,_ i I i_ _ __. __. _;_ i: V: i X A\ - \: Em f:\00".' ~ \ \ \. Figure oi Figure 3-12. Computer olution or Two-Dimensional eat Problem. Figure 3~~~~-12 ontutrSlto:: r — Tl~-i-Diesoa -!ea'r

~AIR-~~~1' ___UNIVERSITY OF MICHIGAN Page 4 CHAYTER 4 SOLUTION OF TEE WAVE EQUATION 4.1 Equation to be Solved One of the most importaint partial differential equations met in physics or engineering is the wave equation. If we let 0 represent the magnitude of a disturbance in aEny medium in which wave propagation can take place, then we can write the wave equation as 2 2 Q*-l) or in terms of Cartesian coordinates x, y, and z +^ ^| (4-2) 3^ 3" a' dz v bt Here v is the wave velocity in the medium and t is the time variable. Equation (4-1) or (4-2) must of course be subject to spacial boundary conditions and initial time conditions. The spacial derivatives of the wave equation have exactly the same form as the heat equation, but the time derivative is second order instead of first order. The dif:ference techniques for converting the partial differential equation to a system of ordinary differential equations in time are also practically identical. Because of this similarity it was decided in this report not to solve any problems involving the application of the electronic differential analyzer to the wave equation, but rather to limit the discussion to a theoretical one showing what the analyzer should be able to do. In order to have a specific problem involving the wave equation to use as an example, the classical problem of the vibrating string will be treated.

.-... ^ a.. | ENGINEERING RESEARCH INSTITUTE Page Page...........A" UNIVERSITY OF MICHIGAN 4 4.2 The Wave Equation for the Stretched String Consider the string showl in Figure (4-1). Let the string be fastened at x = 0 and;> = L, and assumre that the transverse displacement of the string is 0(x,t). If T(x) is the longitudinal tension in the string at (x) = MASS/UNIT LENGTH TENSION c T'/ (xt)= DISPLACE- / /./_ | MENT /. x =L Figure 4-1. Stretched String any distance x, and X/ (x) is the mass per unit length of the string at any x, then the equation of motion for the string becomes (T a a x L 2 (4-3) at with boundary conditions 0(0,t) = 0(L,t) = 0. (4-4) The left side of equation (4-3) represents tie transverse force on a unit length of the string due to the curvature in the transverse displacement; the right-hand side represents the inertial force on a unit length of the string. If the tension T in the string is a constant independent of x, then equation (4-3) can be written a 2 =7; 2

AIR- 1 ENGINEERING RESEARCH INSTITUTE._ A _..R..l E. UNIVERSITY OF MICHIGAN Page which is seen to be the wave equation (4-2) in one space dimension. Thus, the velocity of wave propagation in the string will be given by v - j (4-6) and in general may be a function of x. When the wave velocity is a function of spacial coordinates, equation (4-5) is much more difficult to solve by separation of variables, and hence it is of considerable interest to have a direct method of solution by means of a computer. The initial conditions for our string problem can be written in general as 0(x,O) = f(x) (4-7) (x,O) = g(x) t 4.3 Derivation of the Difference Equation for the Stretched String We now proceed to break the string up into N cells a distance x apart. The difference equation for 0n' the displacement at the nth station, is given by 2 dt n+l -n n-l V dt n Here v = T n, where vn is the wave velocity at the nth station and/ n is the mass per unit length at the nth station. The boundary conditions given in equation (4-4) are imposed by letting 0o = 0N = O. The complete set of difference equations then becomes V2 l2: - 2 - 20

AIR-1 ENGINEERING RESEARCH INSTITUTE ae 46 A UNIVERSITY OF MICHIGAN ( x-.. 2 3 4 0 3 02 3 (4-9) (i )2 N-3 = ON-2 - 2N-3 + 0N-4 (N-3 2 N-2 0N-1 - 20N-2 + N-3 ( _22 N-N-1 N-22 (Nx)2 N- 2ON1 + 2 The equations are subject to initial conditions 01(0), 02(0),... N1(0) and 01(O), 02(0), e. N1(). For some types of problems involving the wave equation, 0/S x may be specified at the boundary instead of 0. In this case the boundary occurs half-way between two stat~ions, and the boundary condition is imposed by equating the 0's on either side of the boundary (See Section 3.3). 4.4 Computer Circuit for Solving the Stretched String Although the set of difference equations (4-9) was derived for the stretched string of Figure 4-1, it clearly represents the wave equation (4-1) in general when the propagation is in one direction. The electronic differential analyzer circuit for solving equation (4-9) is shown in Figure 4-2. 4.5 Solution of the String Equation by Separation of Variables Once again it will prove instructive to solve for the normal modes of vibration of our stretched string by separating time and space variables. By comparing the normal mode shapes and frequencies obtained by this method with those gotten from the difference equations, we can evaluate critically the accuracy of the difference method as a function of the number of cells.

ENGINEERING RESEARCH INSTITUTE Page -Iirt-~1 I_ _ UNIVERSITY OF MICHIGAN I I (Ax)2 i(Ax)d-IL05 (Ax) I1Sx2 2 ~o. AX).5 Ax) 2 2 W W 62 M ^ ^ V2 2 (A) (x)2.. Fvu 3 Copue C f o 0.5 2(I I IN3 N-3 -N-3\ VN -3 (Ax)2 -AAAArI I i I -#N-2 N-2\ i V-2 (Ax)j^^ i r^ /S,^ ^ r^ (Ax)2 I Fiue42 optrCrutfrLovn h aelqa Io 2F~~~~~~~Iu I I:1-1 I% I VN-$r k. ~~iue ici o ovn ie'ae~u o.

ENGINEERING RESEARCH INSTITUTEI AIR-'l UNIVERSITY OF MICHIGAN pp 4 As in Section 3.6a, we begin by assuming that 0(xst) = X(x) Y(t). (4-10) Although it is not a necessary simplification, we make the further assumption that the string is uniform, i.e., / = constant. The wave equation (4-5) then becomes tt It X ~2 Y= X =-O2 = constant (4-11) from which It 2 x (x) + o2 X(x) = 0 (4-12) and T 2 Y"(t) + o 2 Y(t) = 0. (4-13) /' The general solution to the second of these equations is Y(t) = A cos / ( t + B sin / t. (4-14) The boundary conditions of equation (4-4) can be met only if X(O) = X(L) = 0. The latter condition limits solutions of equation (4-12) to the following: X(x) = sin o nx (4-15) where ~n L= n = 1, 2, 3,... (4-16) The general solution to the problem is then 00 0(x,t) > (A cos jnt + Bn sin (nt) s in i L) (4-17) -. n nn cnt L (4-17)n nol where the normal-mode frequency of oscillation 60n is given by Ii- I..I. IIIII-'-;. I~I III _ I I IIL._II.I..III III.,I.III.1 IIl..III -!i

~AITR 1 1 ENGINEERING RESEARCH INSTITUTE.....A....._ UNIVERSITY OF MICHIGAN n - nL,n n = 1, 2, 3... (4-18) The constants Al, A2,..., B1, B2,... are evaluated from the initial conditions (4-7) by the method used in Section 3.6a. Thus, the normal modes are sinesoidal in shape and the frequency of the nth mode is n times the frequency of the fundamental, which is (1T/L) /. We will now solve for the normal modes of the difference equation and compare them with the true normal modes given above. 4.6 Difference Equation Solution of the Stretched String We have already derived the difference equations obtained when the string is divided into N cells [see equation (4-9)] If our string is unifor so that v = /T//J = constant, then we can always select our units of time and distance so that Ax = 1 and v = 1. The length L of the string then becomes L = NAx = N, and we have from equation (4-8) for the ith cell 0i = ji+l - 20. + 0i- (4-19) Let us assume that the normal modes of vibration of the difference equations have the same shape as those for the continuous string. Then for the nth rccc u. = a sin nti sin nt (4-20) and nT i n T + n — i nt ui+ = a(sin N cos c si t (4-21) Substituting equations (4-20) and (4-21) into (4-19) we have 2 = 2(1 - cos n) (4-22) as the expression for the normal-mode frequency ) n. From equation (4-20) it is apparent that uo = uN = 0, so that our boundary conditions are Riet. Equation (4-22) can be expanded in a power series giving

l AIR~ 1 | ENGINEERING RESEARCH INSTITUTE AIR-1 Page 50 UNIVERSITY OF MICHIGAN _ Xn = N()2 l1 12 n ( )2 J (4-23) In the limit of infinitely large N the An above approaches the Wn given in equation (4-18), since here L = N and T = 1. For a finite N the normalmode frequencies from the difference equations are evidently lower than those for the continuous string, A plot of this frequency deviation as a function of the number of cells for the first seven modes is shown in Figure 4-3. 4.7 Summary of Investigation of Difference Techniques Appliedt the Wave Equation Xi ~ The equation of the stretched string has been solved both by separation of variables and by replacing spacial derivatives with finite differences. Comparison of the normal-mode shapes in both cases shows exact agreement. Comparison of normal-mode frequencies shows the difference equation frequencies are somewhat lower than the true frequencies, particularly for higher modes or if fewer cells are used (see Figure 4-3). Since the string equation is the wave equation, the above remarks apply to the solution of the wave equation in general by difference techniques. Thus far we have considered partial differential equations with boundary conditions occurring a finite distance apart. It seems evident that our difference techniques as used here are limited to this type of equation, Thus, it would not seem possible to solve problems in semi-infinite or infinite media unless one can let the time variable in the computer represent the spacial variable which goes to infinityo It should be straightforward to solve problems having spacial coordinate systems other than Cartesian, eogo, cylindrical, spherical, etc. For the appropriate geometries this would undoubtedly require many less cells to realize a desirable accuracy,

ENGINEERING RESEARCH INSTITUTE AIR-1 UNIVERSITY OF MICHIGAN Page 51 NUMBER OF CELLS 2 4 6 8 10 12 14 16 18 20 22 24 LU. z ___ z z b I- I d te C/ /I/ I0 C.)II. w 12 Figure 4-3. Comparison of Norrnal-liode Freqiuencies of thie Difference Equations and the Continuous Jtring.

ENGINEERING RESEARCH INSTITUTE A_^.-I _ |I _ _UNIVERSITY OF MICHIGAN Page CHAPTIJI 5 5.1.juation to be Solved It woula see]i.. of particular engineering interest to investigate the usefulness of' t-he diffierence technique in solving the problein of foleural vibration of beams. lork of this type nas been done on the Caltech analogue 4 computer. In this chapter we will investiate tne theoretical possibilities of the difference technique and ipresent the results of solutions by the electronic differenltial analyzer at tihe University:jf i.ichiaian. ---- ^' 5^' -aw-alw YzDEFLECTI0N'-f>ft___ X I X:U X:L Figure'-1. iibratin, Hea-I. Consider the vibrating beam shown in figure 5-1. If we liirit ourselves to the transverse deflection shoiwn anc. as-eie that the flexual planes re!ain parallel, then for smrall deflections thle equation of motion is given by / 2 2 (x) 9t2 + 2 LI(x)' x2 = 0 (5-1) /t 9x x w;here x = horizontal distance from the left ecnd of the beam t = time y(x,t) = transverse deflection of tile beam at any instant p (x) = mass per unit length of beaimi, at x

ENGINEERING RESEARCH INSTITUTE AIR-1 Page 53 l~ ___AUlnl1 _____UNIVERSITY OF MICHIGAN I(x) = area moment of inertia, at x E(x) = modulus of elasticity, at x EI(x) = flexural rigidity at x. In the derivation of equation (5-1) the effects of rotary inertia and deflection due to transverse shear are neglected. This is valid when cross-sectional dimensions are small compared with the length. We observe that the bending maoment M(x,t) is given by 2 M(x,t) = EI(x) ay(x t) (5-2) 2x2 whereas the shear is given by V(x,t) =- h(x,t). (5-3) ax Equation (5-1) is of course subject to both boundary and initial conditions. The boundary conditions depend on the type of end conditions. Various end fastenings and the appropriate boundary condition at the end are summarized in the following table: Boundary Conditions End of Beam Boundary Conditions Free M = 0, V = 0 Hinged y = O, M = 0 Built-in y = O, ay/ ax = O The two types of beams of most general engineering interest are the free-free beam (both ends free) and the cantilever beam (one end built-in, the other free). We will also consider the hinged-hinged beam because it is the easiest to analyze theoretically and will give us a good idea of how the other beams will behave when difference techniques are used. 5.2 Derivation of the Difference Equation for the Vibrating Beam Once again we will convert the partial differential equation (5-1) for the vibrating beam into a set of ordinary differential equation by using the difference technique. Thus, distance along the beam is broken into N seg-, U J_ - -. ill'I i i ~ I I i i I i! Ij I!iiI

ENGINEERING RESEARCH INSTITUTE Pae AIR-1 Page UNIVERSITY OF MICHIGAN ments of width Ax; the displacement yn at the nth station will then be a function of time only. Following the method of Section 1.2, we have from (5-1) as the equation of motion of the nth cell d2 (ax)?n + Ml - 2M + M = (5-4) n dt2 n+l n h-l where El. M (yi+l - 2yi+ y ) (5-5) 3i (AX) 2 We also note that M - M n Mn-l V-1/2 Ax n- (5-6) and — X~_ - -Sn' Y1nn n-i (57) 1/2 Ax' Before writing down the complete set of difference equations for N cells, it is necessary to consider the boundary conditions. 5.3 Representation of Boundary Conditions in the Difference Equations Assume we have an N cell beam and wish to impose the boundary conditions associated with a particular end fastening, e.g., a free end at the right-hand extremity of the beam. This means that both the shear V and bending movement M must vanish at the beam end. Let us assume, then, that the end occurs at N+l/2 and that VN+l/2 0. From equation (5-6) this implies that N = + = O. But from equation (5-4) this means that 2 (ax) / N 2 MN-1 d. -.., - 1 2 =.'fLN-l dtI

ENGINEERING RESEARCH INSTITUTE _____AIR-1 __ I UNIVERSITY OF MICHIGAN Pe The remainder of the equations are similar to equation (5-4) until.the lefthand boundary is reached, at which point the difference equations a-,ain depend on the type of end fastening. Following the same line of reasoning as above, one obtains the following set of conditions for the difference ecaations for various end fastenings of an N' cell beam: End h'.ere End Occurs Condition Free N+l/2 Il = l"l = 0 Hiinged I 1 =y = 0 Built-In 1N+l/2 y1N y+ = 0 The actual way in wvhich these conditions modify the difference equations is best seen by considering a specific type beam, as in the next section. 5.4 Computer Circuit for holvinjg the Cantilever 3eam by Diffcrence Techniques Since it involves both a free eiind and a built-in end, the cantilever beam shown in Figure 5-2 seeias the best choice'or a specific example. The / Fige 2 3 4 5 6 7 8 Figure 5-2. Cantil evr deam.

AIR-i l ENGINEERING RESEARCH INSTITUTE 56 UNIVERSITY OF MICHIGAN Page left-hand end of this beam occurs at station 1/2, while the right-hand end occurs at station N+1/2. From equations (5-4), (5-5), and (5-7) along with the boundary conditions of Section 5.3 we obtain the following set of difference equations: d2 (Ax)2 2 2 + M3 - 2 + =,2d 2+ - dt2 (Ax)2f3 -2 + M - 2M3 + 2 = ~ * * * * 0 (5-8) (2ax, — t -.-. (Ax) N-2 + MN- 2 + MN-3 = ~ d 2 ( NdX2-1 dt2 - P (A x) PN 2 + MN1 = 0 where EI1 M1 = - 2 Y2 2 (Ax) 3 2 )2 EI.......... ~(y,,i i. i i-.)

ENGINEERING RESEARCH INSTITUTE P AIR-1 Page 57 UNIVERSITY OF MICHIGAN (5-9) EI N2 MN-2 = x)2 (YN-1 2N-2 + YN3) EI N-1 MN-1 (= 2 (YN - 2N-1 + N-2) Notice that even though the left-hand end of the beam occurs at station 1/2, the displacement yl at station 1 is held fixed at zero. The computer circuit for solving equations (5-8) and (5-9) is shown in Figure 5-3. Initial conditions on yn and n must of course be specified in an actual problem. 5.5 Theoretical Solution of the Difference Equations for Vibrating Beams (a) Hinged-hinged beam. In order to check the accuracy of the difference method for beams, we will now solve for the normal modes of vibration of a hinged-hinged beam by separation of variables. The resulting mode shapes and frequencies will then be compared with the difference equation solution of the normal modes for a hinged-hinged beam. From equations (5-1) we have for a uniform beam f t2 + 4 = (5-9) with end conditions 2 2 y(0,t) = y(L,t) = (O,t) = - (L,t) =. (5-10) We separate variables by assuming that y(x,t) = X(x) T(t). (5-11) iiiii~~~~~~a i ii i i i_ ii ii iiiiiIiIiiii

~AIR-1 ENGINEERING RESEARCH INSTITUTE A_____.1 __ IUNIVERSITY OF MICHIGAN Page 58 EII/h1 4'M I~~~ EIA/~dxf A AAA 05.5.5 O(A)X2p5.5 73~~~~N -M4 /\ — ||J_ _111 /^~~~EI3/(Ax)' / ^ / \ A-VW~r-3. I ya by the Difference kethod._34h / 5 \ i^ (L^^N- r^ / 5 NwMr Figur 5-3. Goiue ici o ovn h atlv B a byteDfeec'ehd

~AIR-1 ~ ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN Substituting equation (5-11) into (5-9), we have It IV Tl l X - 2 = constant (5-12) EI T X from which tt 2 EI T + o( - T =0 (5-13) and X + (< =Q 0 (5-14) The solution to the first of these equations is T(t) = A coscK t + B sin</ t. (5-15) The only way the boundary condition (5-10) can be met is if X(0) = X(L) X (0) = X (L) = 0. Hence, equation (5-14) has the solution X(x) = sin /x = sin nrx = 1 2 3, (5-16) sin n i~, 2, 3, (516 from which v( can have only discrete values o( given by - ~ 2 n L (n = 1, 2, 3,... (5-17) Thus, from equations (5-15) and (5-17) the normal mode frequencies W & are 2 2 ElI~-~ ~~~"'2T2 T T (5-18) The shape of the modes is seen from equation (5-16) to be sinesoidal. Let us now proceed to solve for the normal modes of the uniform hinged-hinged beam when the spacial derivatives are replaced by finite differences. If we select our unit of x so that Ax = 1 and hence L = NAx = N, the equation of motion for the ith cell becomes from (5-4) and (5-5)

ENGINEERING RESEARCH INSTITUTE 60., At"' - \ ___,UNIVERSITY OF MICHIGAN d y. dt2 yi+2 4Yi+1 6Yi- 4Y-1+ Yi-2= 0. (5-19) dt We now assume that the normal modes of vibration of the difference equation are the same as those for the continuous beam, namely, sinesoidal. Then for the nth mode yi = a sin N sin A nt (5-20) where n is the normal mode frequency for the nth mode. It is evident that this solution satisfies the boundary conditions that y = M = 0 and YN = TM = 0. Substituting equation (5-20) into (5-19) we find that =n[3 - 4 Cos 2n+ ] (5-21) n\ P N N(5-2 is the expression for Xn. After expanding the cosine functions of (5-21) in a power series, it follows that — 2 EI (n T)4 1 )L2 + ]. (5-22) n QN4 6 N Since here N = L, comparison of equation (5-22) with equation (5-18) shows tkiat as the number of cells N becomes very large, the normal-mode frequency \n approaches the value on for the continuous beam. For finite N the difference equation frequency An is somewhat lower than the true frequency Wn" A plot of the percentage deviation in normal-mode frequency versus the number of cells is shown in Figure 5-4 for the first five modes. (b) Free-free beamui One can solve for the normal modes of vibration of a uniform free-free beam by separating variables, just as we did for the hinged-hinged beam. However, in order to find the normal-mode frequencies w it is necessary to solve for the roots of a transcendental equation. If we define a dimensionless frequency parameter /3 n as

ENGINEERING RESEARCH INSTITUTE Page 61 AHL~~R-1 ~ UNIVERSITY OF MICHIGAN II >* NUMBER Of CELLS 0 4 1 12 14 16 18 24 z,-loM 4 -- - _i /i /r! 1 Igue.. o z -I - - HINGED-HINGED BEAM Figure 5-4. Frequency Deviation for a Hinged-Hinged Beam.

- ENGINEERING RESEARCH INSTITUTE AMR-1 Page 62 -AIR-i^ E__ UNIVERSITY OF MICHIGAN =, bn (5-23) then the values of n for the first five modes of a free-free beam are: 1 2 3 4 5 22.4 61.6 121.0 199.8 298.6 It is not so straightforward to solve for the normal modes of the difference equation representation of the uniform, free-free beam. If we choose our units of x so that A x = 1 and our units of time so that ^ =1, then from equations (5-4) and (5-5) we have for the ith cell i+ Yi+2 - 4Yi+l + 6y - 4yi-1+ Yi-2 (5-24) After we apply the boundary conditions we are left with a set of N ordinary differential equations. Contrary to the case of the hinged-hinged beam, it is not easy to solve for the normal modes here in general. Rather, we have to choose a particular number of cells N, write down the specific difference equations, and solve for the normal modes by a direct but tedious process. Suppose, for example, we decide to solve for the normal modes for 8 cells. If we choose the ends of the beam at the 1/2 station and 8-1/2 station, then from equation (5-24) and the boundary conditions in Section 5.7 for a free-free beam we have the following set of equations: 1 + Y 1 - 2Y2 + Y3 = 0 Y2 - 2yl + 5Y2 - y3 + Y4 = 0 (5-25) + Y3 - 4YL + 6y5 - 4y6 + Y7 = o ~~~~~~~~iiiiiiiiiiiiiiiiiiiiii-iii iiiiiiiiiiI i i...

|R ENGINEERING RESEARCH INSTITUTE 63 _ AIR-i ENGINEUNIVERSITY OF MICHIGAN.. 6 Y6 + 4 - 45 + 6Y6 - 4Y7 + Y8 = 0 Y7 + Y5 - 4Y6 + 5Y7 - 2yg = 0 Y+ 7 + Y 2 + = 0 In order to solve for the normal modes of vibration we assume that at the ith station, yi varies with time as sin t, where, is the frequency of the oscillation. Making this substitution in equations (5-25) gives us (1- X2)y - 2y2 + y3 = 0 -2Y1 + (5- 2) 2 - 4Y3 + y4 = 0 Y - 4y2 + (6 - )y3 - 44 + y5 = 22 5~~~ (5-26) Y2 - 4y3 + (6 - 2 )y4 - 4Y5 + Y6 ~ (5-26) y3 4y4 + (6 - X2)y - 4y6 + y7 =l y4 4y5 + (6 - ) )y6 4Y7 + y 0 Y5 - 4Y6 + (5 - 2)y7 - 2yg = 0 Y6 - 27 + (1- )yg = 0 The only nontrivial solution of these 8 simultaneous algebraic equations is the one for which the determinant of the coefficients vanishes; i.e., we have to eliminate the y's and find an equation in X. When the tedious algebra is carried out, one is left with the following equation in )2 -336 + 3312 X2 - 5140X4 + 2432 X6 - 456 \8 + 36 X10 - X12 = O. (5-27) The roots of the polynomial then represent the normal-mode frequencies. The first three values of Xn obtained from equation (5-27) are

ENGINEERING RESEARCH INSTITUTE Pag _ AIR-1 __ ENGINE__ UNIVERSITY OF MICHIGAN P64 = 0.352, X2 = 0.940 3= 1.74 From equation (5-23) we must evidently multiply A nby 64 to obtain n for the 8 cell free-free beam, since now L = N = 8, and EI/p = 1. A summary of the dimensionless normal-mode frequencies 3n for the continuous beam and the cellular beam is shown below. Free-Free BeamI Mode (continuous beam) (8 cells) % deviation 1 22.37 22,55 + 0.8 2 61.7 60.1 - 2.5 3 121.0 111.3 8.0 Comparison of these results with those for the hinged-hinged beam with 8 cells (see Figure 5-4) indicate that the difference method gives somewhat more accurate frequencies for the normal modes of a free-free beam, Tn order to solve for the mode shapes, it is necessary to substitute the roots Xn back into equations (5-26), In Figure 5-5 a comparison of the 8-cell mode shapes with the mode shapes for the continuous beam is made. ~- ~ (c) Cantilever beam. When the space and time variables are separated in the equation for a cantilever beam and the normal-mode frequencies are determined, the following values of the dimensionless frequency parameter /3n defined in equation (5-23) are obtained: 1 2 3 45 3.516 22o03 61,7 121*0 199.8 If the units are selected so that Ax = 1 and ft/EI = 1, equation (5-24) represents the difference equation for yi when the cantilever beam is broken into cells. Again let us consider 8 cells.'When the boundary conditions outlined in Section 5o3 are applied (left end built in, right end free) the following set of seven simultaneous difference equations result:

ENGINEERING RESEARCH INSTITUTE Page 65 f~Ai~~~IR-1 __UNIVERSITY OF MICHIGAN Ie --------— 11-.8 --— FIRST MODEI-.6 z. _|\ STATIONS ALONG THE BEAM a~~~~~~~.E ~I I 1.0 ------- CELLULAR BEAM 1d *C nV — eaL.4I. - STATIONS ALONG/THE BEAM y, 1 \2 4 5 6 7 8 W-.4.6 Figure 5-5. Comrparison of hode shapes for 8-Cell Free-Free and Continuous Beam.

AIR1 I ENGINEERING RESEARCH INSTITUTE Page 66 n~,-l- [,,,__ UNIVERSITY OF MICHIGAN Y2 + 6Y2 - 43 + Y4 = o 3 - 4y + 6y3 - 4y4 + y5= 0 Y4 + Y2 ~ 4y3 + 6y4 - 4y5 + y6 = 0 (5-28) Y5 + y3 - 44 + 6y5 - 4Y6 + y7 0 Y6 + Y4 - 45 + 66 - 4Y-7 + Y8 = Y7 + 5- 4Y6 + 5y7 - 2y8 = 0 Y8 + Y6- 2y7 + 8 = 0 As before, we assume yi varies with time as sin X t. Equations (5-28) are then reduced to 7 simultaneous algebraic equations. Eliminating the y's gives us 1 - 3362+ 33124- 51406+ 2432X - 456 X10+ 36 >12-. (5-39) 2 The roots of the above polynomial in > are the normal mode frequencies. The first four values of X obtained from equation (5-29) are n 1 0.0554, A 2 =. 347,, )3 0.,940 = 1.77. For our 8-cell cantilever beam Ax = 1 and L = N = 8. Also, EI/p = 1 and hence the dimensionless normal-mode frequency / nof equation (5-23) is obtained by multiplying X by 64. In the following table n for the 8-cell beam is compared with /n for the continuous beam. Cantilever Beam Mode (continuous beam) / (8 cells) % deviation 1 3,516 3.545 + 0.8 2 22.03 22.2 +.0,8 3 61.7 60.1 - 2.5 4 121,0 111.3 - 8,0 L1, 1 i 1 1[i ~ ~I l[ [I [II l........ ~'..... " "!-'-1

AIR-1 I ENGINEERING RESEARCH INSTITUTE Page 67 __ AIR-l"1 __ _UNIVERSITY OF MICHIGAN Evidently an 8-cell uniform cantilever beam (actually requiring only 22 operational amplifiers) gives tolerable normal-mode frequencies for the first four modes. For many engineering problems this would be entirely adequate. Difference-equation representation of a nonuniform cantilever beam would presumably give the same order of accuracy. The mode shapes are calculated in the same manner used for the freefree beam. In Figure 5-6 these mode shapes are compared with those for the continuous beam. Agreement seems to be entirely satisfactory. 5.6 Computer Solution for a 5-Cell Cantilever Beam By choosing our units of x so that Ax = 1 and our units of time so that /EI = 1, then from equations (5-4) and (5-5) the equation for the ith cell becomes i 4+ i- + Yi2 2i+l i i + (5-24) For the 5-cell beam under consideration here let us assume that the clamped end occurs at the 1/2 station and the free end at the 5-1/2 station. Then after applying the boundary conditions as in Section 5.5c we are left with four differential difference equations: Y2 + 6Y2 - 43 + 4 = 0 3 - Y2 + 6y3 - 4y + y5 = 0 +ee00~~~ y2 43 =0(5-30) 4 + Y2 - 4Y3 + 5y4 - 2y5 = 5 + Y3 - 2y + y5 = We note that the length of our beam is L = N Ax = 5 for 5 cells. From equation (5-23) and the table for 3n in Section 5c the period of the fundamental mode for an N cell cantilever beam is T2 2 2 2 N2 T - = ~ -2T...... (5-31) EI n~

,AIR- I ENGINEERING RESEARCH INSTITUTE Page 6 ___IAIR-1 __ _UNIVERSITY OF MICHIGAN 1.0. 0.k FIRST MODE 0.6 - --- 0.4 - B - __"___ _ ____SN TH B CLL 0. 0.0 Z 1.01 I 2 3 4 5 6 7 8 -0.2 - STATIONS ALONG THE BEAM.0.4 — I__ ____ — - | SCONTINUD MODEUS.0.6- BEAM 0.48-...0o.4.I I IEI U Figure 5-6. Comparison of Iode Shapes for 8-Cell Cantilever Beam and Continuous Beam.. 0.2 - 0.8. and Continuous B3eam.

AIR-| I ENGINEERING RESEARCH INSTITUTE P 69 AJIRI UNIVERSITY OF MICHIGAN Pa since EI/p = 1 for our beam. Thus, when N = 5 we find that T1 = 44.7 seconds, which is somewhat inconvenient for the electronic differential analyzer. In order to shorten the normal-mode periods by a factor of five the input resistor to each of the integrating amplifiers is reduced by a factor of five. Then one unit of computer time actually is the equivalent of 0.2 seconds. The computer circuit for the 5-cell cantilever beam is shown in Figure 5.7. It was decided to examine the response of the beam to a step force applied at station 3, i.e., at the middle of the beam. In order to do this, a voltage V(t) was summed with the various feedbacks of station 3 (see Figure 5-7). A constant voltage was first applied and resistors were inserted across several of the integrating condensers to damp out the resulting beam oscillations. The damping resistors were then removed and finally the constant voltage input was released. This corresponds to suddenly removing a force at station 3. The resulting computer response at each of the four stations is shown in Figure 5-8. Note that the second harmonic is clearly visible at stations 2 and 3. Higher harmonics (There are four in all) are visible in Figure 5-9, where the acceleration y2 (near the clamped end of the beam) is shown as a function of time. From Figure 5-8 the periods of the fundamental and second harmonic are 8.90 and 1.42 seconds, corresponding to 44.5 and 7.1 units of computer time. From equation (5-31) the value of 3n for n = 1 and n = 2 can be calculated. In the following table these values of ~1 and 32 are compared with those for a continuous contilever beam and those calculated for a 5-cell beam using the method of Section 5.5c. Mode (continous beam) (computer 5 cells) (theoretical, 5 cells) 1 3.516 3.54 3.59 2 22.03 22.1 22.3 The theoretical values of A for a 5-cell cantilever beam were calculated by assuming yi varies in time as sin X t. Equatibns (5-30) then become algebraic and are satisfied only for those values of ) satisfying the equation \/ 1 - 502 + 75X4 - 18X6 + = 0 (5-32) 2 By multiplying the roots of this polynomial by N = 25 as in Section 5.5c, one obtains the values of frequency parameter /3.

AIR~1 |- ENGINEERING RESEARCH INSTITUTE Page _AIR-i "_ __UNIVERSITY OF MICHIGAN 70 1/6 I y2 1/4 0.2 r 0.2 r 0- 3 l 3I Y22 Y4 1/4 1 11 Y2 -y3 AAAA)r- - 1/4Y3 t Y4 -Y5 4y3 rlm l[___ 1/5 y54 1/2 - — Wr-WW'-Y1 — -I-' Y5 = I/ ~**v-vvvrvv V(t) -y5 FEEDBACK CONNECTIONS OMITTED FOR CLARITY Figure 5-7. Cormputer Circuit for 5-Cell Cantilever Beam.

ENGINEERING RESEARCH INSTITUTE Page 7 AIR-1 UNIVERSITY OF MICHIGAN CHART NO. BL: \ T-/ —':lL-L-' CHll NO' B L -2 THE BRUSH DEPN V l l -_ —-- "1 1: —I I /- II I1 -f11111 -- -— 1 —-1 - 1 - -. —--- II i i i itali A ~'~ 1:4: ~IIVi —\~::~-~ - CHART NO. BL-924 THE BRUSH DEVELOPMENT (._Figure~~~~~~~. Figure 5-8. Computer Solution for 5-Cell Cantilever Seam.

AIR-1 ENGINEERING RESEARCH INSTITUTE Pae UNIVERSITY OF MICHIGAN 7 / THE BRUSH DEVELOPMENT CO,. PRINTED IN US.A. It should be remarked that the input resistors to each summing amplifier were matched to 0.1 percent. This is necessary to assure that accurate differences are obtained. This point is analyzed in the next section 5.7 Accuracy Requiremerid for the Difference Technique (a) Component accuracy. One of the fundamental difficulties encountered when continuous derivatives are replaced by finite differences is that the smaller the interval used, the more accuracy is required in taking the difference. This is particularly true when the order of the derivative is high, such as fourth order. In the case of the vibrating beam, it is evident that the greater the number of cells for any half-wave length of a normal mode, the smaller the differences become and the more critical the accuracy requirements become for the summing resistors. For example, consider the 5-cell cantilever beam of

UNIVERSITY OF MICHIGAN Page the previous section. Suppose all the summing resistors in Figure 5-7 are perfectly accurate except the y3 input resistor to the bottom bank of amplifiers. Let us assume that it is 1 percent high, i.e., 1.01 instead of 1 megohm. Then, the difference equations are exactly the same as those given in (5-30) except that the last equation is now Y5 + l.Oly3 - 2y4 + Y5 = 0 The characteristic polynomial for this perturbed system is 1.11 - 49.8?772 7+.99,4 - 18xi + i78 = 0 (5-33) compared with 1 - 50.2 754 _-. I 6 + = 0 (5-32) for the unperturbed system. The value of l obtained from the first root of equation (5-33) is 3.79, which is 5.6 percent higher than the value of 3.59 calculated for the unperturbed system. Thus, a 1 percent error in 1 of the 14 summing resistors of Figure 5-7 results in a 5.6 percent error in the frequency of the fundamental mode. The error in the higher modes will be very much less, of course. This simple example is graphic demonstration of the necessity in maintaining high accuracy in the summing resistors. The main requirement is that the relative values of any set of input resistors to one amplifier be accurate. The absolute values are relatively unimportant. Also, the tolerance in the resistors representing EI andf in Figure 5-3 is not critical. For the summing resistors requiring high tolerance it is planned to use wire-wound resistors matched to the order of 0.01 percent. (b) Effect of small voltage transients. The deflection of a beam due to an applied force is proportional to the fourth power of the length of the beam. The deflection of the beam at any point is represented by a voltage output of the electronic differential analyzer, while an applied force is equivalent to a voltage input to a summing amplifier in Figure 5-7. This latter voltage may be deliberately introduced through an input resistor or may be inadvertently introduced through a sudden shift in the balance of a summing amplifier. As a result, transient oscilla

ENGINEERING RESEARCH INSTITUTE Page 7 UNIVERSITY OF MICHIGAN tions in the whole computer circuit can be caused by fluctuations in power supply voltage. The sensitivity of the network to small voltage inputs will increase as the fourth power of the number of cells. Thus, it becomes important to have an extremely well regulated power supply. Drift stabilized amplifiers would not help a transient condition but would serve to maintain the static deflections to a high degree of accuracy. 5.8- Summary of Difference Technique for Vibratin Beams The vibrating-beam equation has been solved with the difference method both theoretically and with the electronic differential analyzer. Results show that normal-mode shapes and frequencies exhibit good agreement with continuous beams providing 3 or more cells per half-wave length of the normal mode are used.!ji If a large number of cells is used, it is necessary to realize extreme precision in the summing resistors in order to be able to take small differences accurately. Also, a large number of cells means that transient voltages, such as might be introduced by power supply fluctuations, must be held to a minimum, Cantilever and hinged-hinged beams have a fixed equilibrium positio relative to their surroundings and hence are stable on the electronic differential analyzer. Free-free beams are not supported, however, and will tend to be unstable, since any small voltage unbalance will cause them to rotate and translate as well as vibrate. Damping can easily be included in the beam equation by placing the appropriate resistors across integrating condensors. Any variable force as a function of time can be introduced at any point or points along the beam. The final computer response gives directly the bending moment and displacement as a function of time and distance along the beam. This same difference method can be used to solve beams with both torsional and lateral bending.4 In this case the torsional equation is similar to the wave equation treated in Chapter 4, The proper cross-coupling resistors then tie the two systems together.

ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN Page BIBLIOGRAPHY 1. Hagelbarger, Howe, and Howe, Investigation of the Utility of an Electronic Analog Computer in Engineering Problems, External Memorandum UMM-28 (April 1, 1949), University of Michigan Engineering Research Institute, AF Con. W33-038-ac-14222 (Project MX-794). 2. Howe, Carl E., Further Application of the Electronic Differential Analyzer to the Oscillation of Beams, External Memorandum UMS-47 (June 1, 1950), University of Michigan Engineering Research Institute, AF Con. W33-038-ac-14222 (Project MX-794), 3. Howe, Howe, and Rauch, Application of the Electronic Differential Analyzer to the Oscillation of Beams, Including Shear and Rotary Inertia, External Memorandum UMM-67 (January 1951), University of Michigan Engineering Research Institute. 4. McCann and MacNeal, Beam-Vibration Analysis with the Electric-Analog Computer, Journal of Applied Mechanics, Vol. 17, No. 1, March, 1950. 5. Timoshenko, S.P., Vibration Problems in Engineering, D. Van Nostrand.........................,,,,

UNIVERSITY OF MICHIGAN 3 9015 03026 8612