THE USE OF COMPUTERS IN AERONAUTICAL ENGINEERING EDUCATION Robert M. Howe and Gabriel Isakson Department of Aeronautical and Astronautical Engineering The University of Michigan Janu-ary. 1, 1963 This material is distributed by The Ford Foundation Project on the Use of Computers in Engineering Education at The University of Michigan. This report appears in the library edition of the Final Report of the Project and is also issued as a separate booklet. Similar "Curriculum Reports" for other engineering disciplines are also available on request.

ABSTRACT During the past few years the College of Engineering of The University of Michigan has been exploring the use of electronic computers in its undergraduate curricula. This report describes the experiences and current philosophy of the faculty of the Department of Aeronautical and Astronautical Engineering concerning classroom computer use. Included are a description of the Department's curriculum and computing facilities, and a discussion of the extent to which computers have been used in Departmental course work. Seven computer-oriented example Aeronautical Engineering problems with complete analog or digital computer solutions are also included. This set of problems may be considered as a supplement to the 104 example engineering problems, including several of interest to aeronautical engineers, which have been published previously by the Project on the Use of Computers in Engineering Education. -Hl

Table of Contents Page I. Introduction H3 II. The Curriculum H4 III. Analog Computer Facilities in the Department of Aeronautical and Astronautical Engineering H5 IV. Example Problems H5 105 A Third Order Control System H6 106 A Multiple-Loop Aircraft Control Problem H10 107 Solution of the Linearized Longitudinal Flight Equations H19 108 Solution of the Linearized Lateral Flight Equations H22 109 Two-Dimensional Re-Entry Trajectories for Lifting Vehicles H26 110 Free Vibration Characteristics of a Uniform Cantilever Beam with Elastic Root Restraint H32 111 Re-Entry Flight Path of a Lifting Vehicle H43 Table IH Example Problems H5 -H2

USE OF COMPUTERS IN AERONAUTICAL ENGINEERING EDUCATION I. INTRODUCTION The extensive use of computers, both analog and digital, in the field of aeronautics and astronautics raises the question of the role they should play in the education of engineers for work in this field. The present report describes the experience and current philosophy of the Aeronautical and Astronautical Engineering Department of The University of Michigan in this regard Since 1948, electronic analog computers have been used extensively for instruction in the Department. Originally used for laboratory simulation of control systems in graduate courses in automatic control, analog computers have more recently been used for classroom lecture demonstrations of air-frame dynamics, satellite orbits and re-entry trajectories. In addition, they continue to be used in a number of laboratory courses, both graduate and undergraduate, to simulate linear and non-linear control systems. Thus, analog computers have been used both in lecture demonstrations and in the laboratory. When used in the laboratory there are usually only two students operating each computer. The students learn something of the capabilities of analog computers, develop a "feel" for the dynamic behavior of various systems through analog simulation, and gain some experience at analyzing recorder output traces. In the example analog computer problems which are included in this report, it is assumed that the reader has a knowledge of the working principles of analog computers. The problems presented include a number of control-system and airframe dynamics problems as well as a satellite orbit and re-entry trajectory problem. These examples are suitable either for lecture demonstrations or for individual laboratory experiments. The Department of Aeronautical and Astronautical Engineering has placed little emphasis on the application of digital computers in its undergraduate program. Instruction in digital computer programming has not thus far been made a requirement of the program, although most other engineering curricula at The University of Michigan do have such a requirement. It is felt that such instruction would be most meaningful if provided in conjunction with instruction in numerical analysis. With this in mind, the Department is deferring a decision in this area until a current study of the mathematics curriculum in the College of Engineering is completed. The main reason for the reluctance of the Department of Aeronautical and Astronautical Engineering to emphasize digital computer programming is concerned with the manner in which digital computers are used in the aerospace industry. Typically, a company in that industry will have a large computer facility with a staff of programming specialists who handle most, if not all, of the programming work of the company. Under these circumstances it is not essential for the engineering staff to know programming. If such knowledge is desirable in individual cases, appropriate instruction is usually readily available within the company. These circumstances, and the great pressure to provide more room in the curriculum for the teaching of scientific principles,explain the current policy of the Department. -H>

It is recognized that the advent of greatly simplified programming languages may alter the mode of operation of computing facilities in industry and elsewhere and permit more direct access to the computer by the engineer. This, in itself, would call for a reappraisal of the need for instruction in programming at the undergraduate level. At the graduate level considerable use has been made of the digital computer by students engaged in research or thesis work. It has been found that when the need for such use arises, the student can usually acquire the necessary programming knowledge in a sufficiently short time. The areas in which it is presently felt that the most effective use can be made of digital computers are structural analysis and space flight mechanics. In the former, the advent of the digital computer has stimulated the development of numerical methods for the analysis of complex structures. Much of the structural analysis performed in the aerospace industry at the present time involves the implementation of such methods on digital computers. In the field of space flight mechanics, the digital computer shares with the analog computer the status of an indispensible tool in the solution of trajectory problems. Two digital computer example problems are presented in this report. While they have not been used for classroom instruction, it is felt that they would be suitable for that purpose. One of the problems is in the field of structural dynamics while the other is in the field of space flight mechanics. II. THE CURRICULUM Freshman Year (First Semester): Hours Freshman Year (Second Semester): Hours Mathematics 4 Mathematics 4 Chemistry 4 Chemistry 4 English 4 English 2 Engineering Graphics 3 Engineering Graphics 2 Physics 5 Sophomore Year (First Semester): Sophomore Year (Second Semester): Mathematics 4 Mathematics 4 Physi'cs 5 Chem.-Met. Engineering 3 General Aeronautics 1 English 2 Statics 3 Mechanics of Materials 4 Elective 3 Laboratory in Mechanics of Materials 1 Thermodynamics I 3 Summer Session Elective 3 Circuit Analysis and Electronics 4 Junior Year (First Semester): Junior Year (Second Semester): Mathematics 3 Flight Structures I 3 Aerodynamics I 4 Propulsion I Dynamics 3 **Mechanics of Flight Structural Mechanics 3 Aerodynamics II Elective 3 Elective 3 Senior Year (First Semester): Senior Year (Second Semester): Flight Structures II 4 *Automatic Control Systems 4 Propulsion II 4 Airplane Design I 4 Electives 6 Electives 6 Economics 3 English 2. This course has used the analog computer in the laboratory. **e This course has used the analog computer as a lecture demonstration tool. -H4~

Use of Computers in Aeronautical Engineering Education III. ANALOG COMPUTER FACILITIES IN THE DEPARTMENT OF AERONAUTICAL AND ASTRONAUTICAL ENGINEERING At the present time some 13 small table-top analog computers are available within the Department; eight of these have 6 amplifiers each, while the remaining five have 10 amplifiers each. Plug-in quarter-square multipliers and diode function generators are available for use with the 10-amplifier computers. In addition, analog equipment totaling some 116 amplifiers is available with removable patch boards for solving more elaborate problems. Almost all of the equipment has been designed and built within the Department, although any future equipment will probably be purchased commercially due to the availability of suitable low cost computers from several computer manufacturers. It is interesting to note that the original computers used back in 1948 employed external plug-in components, whereas the more recently constructed computers within the Department all have internally mounted computing components and have been designed for maximum programming simplicity and reliability. It is felt quite strongly that these features are important in order to make certain that the student gets maximum utilization from his time in the laboratory and accomplishes his problem-solving goals. Too often the students (and the instructors, for that matter) can get discouraged when attempting to use unreliable or difficult-to-program computing equipment. IV. EXAMPLE PROBLEMS The remainder of the text of this report is made up of seven example problems. Five of these employ analog computers and are suitable to be used either as classroom demonstrations or laboratory experiments. The two remaining problems involve digital computer solutions to problems in space flight mechanics and structural mechanics. The problem titles are listed in Table IH and may be considered as a continuation of problems 1 through 115 published in previous reports of the Project on the Use of Computers in Engineering Education. TABLE IH Example Problems Number Title Author Page 105 A Third-Order Control System R. M. Howe H6 106 A Multiple-Loop Aircraft Control Problem R. M. Howe H10 107 Solution of the Linearized Longitudinal R. M. Howe H19 Flight Equations 108 Solution of the Linearized Lateral R. M. Howe H22 Flight Equations 109 Two-Dimensional Re-entry Trajectories for R. M. Howe H26 Lifting Vehicles 110 Free Vibration Characteristics of a Uniform G. Isakson H32 Cantilever Beam with Elastic Root Restraint 111 Re-Entry Flight Path of a Lifting Vehicle G. Isakson and H43 J. L. Lemay -H15

Example Problem No. 105 A THIRD ORDER CONTROL SYSTEM by R. M. Howe Problem Statement Consider the feedback control system shown in Figure 1. Represent the system by means of an analog computer circuit, investigate the step response of the system as a function of K and correlate the results with the theoretical ones predicted by the Nyquist and root-locus techniques. System Controlled o0 Controller i On 1 01 6 - y -; —= p(0.2p+l)(0.5p+l Input y lOutput Figure 1. Block Diagram of Third-Order Control System Solution The system being controlled has the transfer operator, Q00~ 1 - = G(p) = p(.2p+l)(0.5p+l) (1) y and hence can be described by the differential equation, 0.1'0 + 0.7 o + = y (2) where, when the feedback loop is closed, y = KE = K(Oi- 0). (3) An analog computer circuit for representing the system is shown in Figure 2. An alternative circuit which synthesizes the transfer operator G(p) directly by cascading three separate transfer operators is shown in Figure 3. Here amplifier 2 has an equivalent feedback resistor -=r^ Q^ 0.1. + 7+~ Figure 2. Computer Circuit for Third-Order Control System -H6

equal to 0.1 megohms times the reciprocal of the pot setting, i.e., 0.2 megohms. Thus the amplifier behaves like a simple first-order system with the time constant RC = 0.2 seconds and static gain of 0.2. Similarly, amplifier 3 represents a simple first-order system with time constant RC = 0.5 seconds and static gain of 0.5. The circuit of Figure 3 bears more of a oneto-one correspondence with the physical system being controlled, assuming it consists of the three cascaded elements shown. 1 0.1 0 1 1 1 0.5 1 1/K 1 1 1 1 1 0~1 1 -G. ^y=K~ v ^K —--- -01- 0 0 01 0. o5 _iq 0.2 p+l 0.5 p+l p Figure 3. Alternative Circuit for Third-Order Control System Let us examine the step-response curves of the control system for various gain-constants, K. The results are shown in Figure 4. Note that for K=7, the system is essentially neutrally stable with a frequency of approximately 3 radians/second. For larger K the system becomes unstable. For K a 1 the response exhibits a reasonable transient. The student must be reminded that for large values of K the step input voltage Gi must be kept small enough to avoid overloading amplifier 1, unless the scaling of both amplifiers 1 and 2 is changed appropriately in Figures 2 or 3. K= 0.4. K 106 0o s o o s 0 K=7 A M'=2,5 o 5 lo 0 5 \o Figure 4. Step Response of Third-Order Control System for Various Gain Constants K Having observed the step-response characteristics, let us correlate the results with theory. First, consider the Nyquist plot of G(p) shown in Figure 5. This plot can either be made analytically by computing G(jw) from Eqn. 1 for various &, or it can be obtained from the computer -H7

A Third Order Control System circuit of Figures 2 or 3 by driving the input y to amplifier 2 from a sine wave generator. By recording y and 90 on two channels, the amplitude ratio and phase shift at each frequency can be measured. Output drift of the open loop system can be eliminated by recording 0.1 0 from amplifier 3 and correcting the reading to g0 (900 additional phase lag, gain factor of 10/o). In practice the student should use both the analytic and experimental technique to obtain and plot G(jc). - p plIcne G plane -5 -\ I \ = 3. t real a~x/s O<c K 7, N=O, Z=0, Stable 5 -4 - 7<K, N=-2,Z=2, Unstable Kc0, N=-l,Z=l, Unstable /0,. / K=7, Neutrally Stable, 3 = 3.2 rad/sec Figure 5. Nyquist Plot for the Third-Order Open-Loop System 1 G(p) =p(0.2p+l)(0.5p+l) From the Nyquist plot it is evident that the system is stable for 0<K<7 because the critical point (- -,0) is not encircled, while for K>7 it is unstable. For K=7 we have neutral stability at 3.2 radians/second. This corroborates our computer response curves of Figure 4. Finally, consider the root-locus plot shown in Figure 6. Here for K=O the closed-loop characteristic roots, A1,A2, and 3' begin at the three zeros,Pl=0, P2=-2, and P3=-5, respectively. X3 moves out to the left on the negative real axis for increasing K, while X and X2 move together along the real axis, joining at K=0.4. For higher K they move out into the complex plane, crossing the imaginary axis at K=7 (neutral stability). Again the results confirm the computer records of Figure 4. For K= 1.06 the root locus diagram predicts a damping ratio of g 0.5 and an undamped natural frequency, n - 1.38, radians/second for the oscillatory part of the transient. The corresponding damped transient period of 2w = 5.2 seconds is Mn W -H8

Example Problem No. 105 evident in Figure 4. The third root, \C = -5.6,gives a damped exponential with a time constant of (5.6)-1 = 0.18 seconds. This part of the transient is not evident in Figure 4 for K=1.06 because it is so small in initial amplitude and decays so fast. 10(p+l) G(p) = p(p+2):(p+5) For K=2= XA1 2 n = 5 rad/sec, = 0.62 t \ 3_~ -1/A = 1.25 sec - 5 j2 K= 0,32 —.2 1 Figure 6. Root Locus Diagram for the Third-Order System With G(p) = p(p Discussion In this example control-system problem the student is given an opportunity to set up a control-system simulation on the analog computer, to take recordings and analyze them, and to compare results with the Nyquist and root-locus techniques. -H9

Example Problem No. 106 A MULTIPLE-LOOP AIRCRAFT CONTROL PROBLEM by R. M. Howe Problem Statement Consider the control of the attitude angle of the aircraft shown in Figure 1. Assume the aircraft remains in symmetrical flight (no bank angle, sideslip angle, or yaw angle) and denote the attitude as measured from the horizontal reference by the symbol G0. The aircraft is controlled by the elevator displacement angled. If the angles G0 and S are kept small, a linear differential equation can be used to describe the airframe system and hence a transfer operator Y4(p) = G ^. Figure 1. Geometry of Aircraft System The elevator deflection anglers,will depend on the torques applied to the elevator. These will include the torque from the control-surface actuator (usually hydraulic) and the aerodynamic forces. The latter will depend in a complicated fashion on the airframe attitude angle G0 through the transfer operator Y5(p) as shown in the block diagram of Figure 2. In order to minimize this effect the elevator is positioned by a closed-loop system with inputsi, the desired elevator displacement angle from the autopilot, and with output,6, the actual elevator angle. Let us write the output of the closed-loop elevator servo in terms of the transfer operators of Figure 2. Thus Y2Y Y = i+Y3 si + 3tY2Y (1) where Ta is the aerodynamic torque acting on the elevator. Inspection of Eqn. 1 shows that if we make the servo amplifier transfer operator Y2 sufficiently large, the Ta term in Eqn. 1 can be made quite small and the aerodynamic feedback loop through Y5(p) can then be neglected. Thnis assumption, valid in practice, simplifies considerably the design of the overall control system. Next let us assume that the open-loop transfer operator of the elevator control system takes the form -H10

= Y2y3 p p ) (2) i.e., behaves as a conventional open-loop second order system with velocity constant Kv and time constant F. This implies that the system has inertia, output or viscous damping, and proportional control. The closed-loop transfer operator Y6(p) then can be written 1 Y6(p) 1 p2+ p~(3) + p + p+ 1 (On2 n ne where e = 2 /KT (4) For actual numerical constants we will assume the following values for the undamped natural frequency,/ne, and damping ratio T of the elevator servo: wn = 25 rad/sec, e = 0.4 (5) S^ ervo Amplifier Elevator ~0 Autopilot and Actuator Surface Airframe 0 +i +Se irmT l1 Y 2 I;~ — - Y 3 4 Actuator Torque + Ta Aerodynamic Moment Y5 Figure 2. Block Diagram of Aircraft Control System Although our second-order system model for the closed-loop elevator servo is never exactly correct in practice, it is close enough to a reasonable representation to make our illustration meaningful. Next consider the airframe transfer operator relating attitude angle 90 and elevator displacement S. To simplify the problem somewhat we will assume the aircraft velocity remains constant. This neglects the phugoid motion, which exerts only a minor influence on the transient behavior of the overall control system. By summing pitching moments acting on the airframe, one can show that the transfer operator takes the following form: -Hll

A Multiple-Loop Aircraft Control Problem o0 k(1sp + 1) (6) - Y3(p) P= (a2p + alp + a0) where the constants Ts, a2, a1, and a0 are normally all positive. However, we will assume in our example that a0 is negative. This makes the aircraft longitudinally unstable; e.g., when it pitches up, the aerodynamic moment is such as to make it tend to pitch up further. When the center of gravity of the aircraft is located behind the aerodynamic center of pressure,,such an instability occurs. Most ballistic missiles are aerodynamically unstable in this manner. Let us assume the following numerical values for the constants in Y3(p). 4(p + 1) O Y3(p) = 2 — =- (7) p(o.3p + p - l) 6 from which the equation of motion for the airframe becomes 0.3 go + 90 - G0 = 46 + 4 S (8) Represent the system by means of an analog computer circuit. Add proportional control and determine the effect of the controller constant,B, on the stability of the system. Examine the system theoretically using the Nyquist and root-locus methods and try to predict the observed results. Recommend additional autopilot compensation which should improve the stability and dynamic performance of the system. Solution Usual straightforward representation of this equation on an electronic differential analyzer would require differentiation of the voltage representing 6, which is undesirable. A'better approach is to consider the transfer operator Y3(p) in Eqn. 7 as broken into two parts. Thus Y3(p) = 3 (p) + p3 (p) (9) where 4 90 Y3(p) 3 2 (10) 0.3p + p - p Then 9o = Y3(p)S and o0 = (Y3(p) + 3(p)) (11) From Eqn. 10, 00 satisfies the equation 0.3 00 + 00 = 0o +z (12) which can be set up directly on the analog computer with both G0 and 00 available as integrator output voltages. By adding G0 and 00, the final output 00 is obtained without requiring the -H12

Example Problem No. 106 differentiation of. This is the technique used in the circuit diagram of Figure 3, which also includes the dynamic representation of the elevator servo in accordance with Eqns. 3 and 5, along with an autopilot with constant gain,B,(i.e., proportional control). In the circuit, amplifier 1 represents the autopilot of Figure 2, while amplifiers 2, 3, and 4 simulate the elevator closedloop servo (the aerodynamic moment T is neglected as explained earlier). Finally, amplifiers 5 through 10 represent the airframe dynamics. The block diagram of the simplified system is shown in Figure 4. Computer step-response recordings for the entire system are shown in Figure 5 for various autopilot gain constants B. Note that the response is stable only for 0.385< B <4.7. For autopilot gains lower or higher than this range the closed-loop system is unstable. This is confirmed by referring to the root-locus diagram of Figure 3.6. For large B care must be exercised not to let the amplitude of the input step-voltage 9i be large enough to saturate amplifier 1 in Figure 3. Reference to Figure 4 and Eqns. 3, 5, and 7 shows that the open-loop transfer operator G(p) for proportional control is given by 1 4(p + 1) G(p) = - Yo(p) = 1 2 1). 2 (13) p(-p2 + 0- p + 1)(0.3p2 + p - 1) 1 0.1 l 1 00.12 301 0. 0. 1 0.25 0.1 0.2 01 11 21 2 Qv 1 -AA1 2 3 4>jyU 4 0 0. [/B i.. ^ 6'/62^ E/25 =-6 0.333 1 1 0.1 0.3G0O0+ -0.391 1 0 0 _1 ~1 1 10 3e0.390+00 ~~-H13Figure 3. Computer Circuit for the Multiple-Loop Aircraft Control Problem, Proportional Control G0 Autopilot Elevator Servo Airframe Figure 4. Block Diagram of Simplified Aircraft Control System -H13

A Multiple-Loop Aircraft Control Problem 3 t0 t o7Yki t V - X B=.35 V.385 a. a 0 5 02,5 t E3 = 4.7 | B=5 Figure 5. Step Response for Autopilot System with Proportional Control for Various Gain Constants, B. The Nyquist plot for G(p) is shown in Figure 7. Since G(p) has a pole at p = 0, the p-plane contour makes a small counter-clockwise semicircle about the origin. The corresponding G-plane contour is a large clockwise semicircle as shown in the figure. Reference to the figure shows that the G plot crosses the real axis at the point -2.60 for 0 = 1.3, and again at the point -0.213 for X = 7.5. There are 4 ranges in autopilot gain factor which give different encirclements of the critical point (-l/B,O). These are summarized below. Range in B N 0 <B <0.385 -1 0.385 <B <4.7 +1 4.7<B -1 B<0 0 -H14

Example Problem No. 106 j25 \ B-4.7 \// \ I,/3 \2 /3 -25 -Zo 15 -o - ~040 3 5 t-o / \ 5.\5 /,a s -j 2 Figure 6. Root Locus for Airplane Autopilot System With Proportional Control Since N = P - Z and we wish to determine Z, the number of roots of 1/B + G lying in the right half p plane, it is important to determine P, the number of poles of G lying within the right half p plane. This can be done by inspecting the zeros of the denominator of Eqn. 13. The zeros of the quadratic term representing the elevator ()2 + 8p + 1 are both in the left plane. The zero at p = 0 represents a pole of G(p) which we avoided in our contour in the p plane. The quadratic 0.3p2 + p - 1 has two zeros, namely 0.80, -4.14. The first of these does lie in the right-half p plane, so that P = 1. This pole of G(p) in the right half p plane accounts for the airframe pitching instability. Thus we are dealing here with an unstable openloop system, as reflected by the fact that P = 1. For a stable closed loop system Z must equal zero, and since N = P - Z, N must equal 1 for stable closed loop operation. But from the Nyquist plot of Figure 7 it is evident that N = 1 only if 0.385<B c4.5. For all ranges of B the following values of Z are obtained: -H15

A Multiple-Loop Aircraft Control Problem WA 0.4 / \p planeA / \'^ \ 3 I \ /1 I X07 / I -~6. 0 - o - -.. ret/l axis For\ 0 385 we have neutral closed -loop stabi y'w h -.o 3 F \ / A- 7.5, // ~ (-o.24,o0) \ / \ / K ^~~~/ \ // // \ i \ \ I,' G plane \ / \ Figure 7. Nyquist Plot for Airplane-Autopilot System Range in B N Z 0 <B <0.385 -1 2 0.385 <B <4.7 +1 0 4.7<B -1 2 B <0 0 1 For B = 0.385 we have neutral closed-loop stability with X = 1.3 rad/sec. For B<0.385 there are 2 characteristic roots of the closed-loop system with positive real part, and the system is unstable in an oscillatory manner. For B = 4.7 we again have neutral stability with MO = 7.5 rad/sec, while for B>4.7, Z = 2 and the system is unstable in an oscillatory manner. For B<0 (i.e., negative) Z = 1 and must represent a positive real characteristic root. The resulting instability is pure exponential. This is an example of a conditionally stable system, i.e., one which goes unstable for too low a gain constant as well as too high a gain constant. Step response curves for the closed-H16

Example Problem No. 106 loop system with different gain constants B were shown in Figure 5. Our determination of the range of gain constant B for closed loop stability using the Nyquist method agrees with the rootlocus plot of Figure 6. Note that the response in the stable region is still quite oscillatory. To improve performance let us modify the autopilot transfer operator Yl(p) from pure proportional control by adding bandwidth-limited error-rate control with error-rate constant C = 1 second. Thus let e Y1(P) (+O.l(p) (14) where the time-constant in the denominator is 0.1 C = 0.1 second. Now the entire open-loop transfer operator G(p), written in factored form, becomes 1 83,300 (P+l)2 G(p B =O ( p) pp(p+10-j22.9) p+10+j22.9) p4.4)(p-0.80) (15) The new root-locus plot is shown in Figure 8 and predicts reasonable closed-loop performance for 1/3 <B<2/3. This is confirmed by the computer step-function response curves of Figure 9. The analog computer circuit used to represent the autopilot transfer operator of Eqn. 14 can be synthesized by rewriting Y6(p) in the following manner: Y6() = 1p = 1 + O.- (16) 1+0.lp l+0.lp The appropriate analog circuit of Figure 10 replaces amplifier 1 in Figure 3. Again, care must be exercised to avoid step inputs which are large enough to cause saturation of amplifier 12 for large gain constants B, unless the whole circuit of Figure 3 is rescaled. \ 3 / B=0.61 9=0.937 \ I/ / 9=0.7 - / \/ / \.\ 47 / \ - \ i-I~72~ Ii / V-..5 / \ \ / /\ \ /3=2.3.B456 \ Figure 8. Root Locus for Compensated Autopilot System -H17

A Multiple-Loop Aircraft Control Problem e..5. /.o B —2 o, lo 15 0 t 10:B 1.72 Figure 9. Step Response of Compensated Autopilot System for Various Gain Constants,B. "00 ~~~IC~~~~~~~~~ DI? 1 1 >1. i 1 - i0.1 12-0.9 0 -.5'I 0 "A!.I..-. —.A.-.j r X B 11 Figure 10. Analog Circuit for Representing the Autopilot Transfer Operator Y6(p) l+. p = i-o 1+0.Ip 8-H18-00 -H18

Example Problem No. 107 SOLUTION OF THE LINEARIZED LONGITUDINAL FLIGHT EQUATIONS by R. M. Howe Problem Statement Let us consider the linearized longitudinal flight equations of an aircraft. They can be written in terms of the dependent variables, velocity perturbation, vp, pitch-rate, q, and angle of attack perturbation, a. These equations can be obtained by considering perturbations from body axes which are displaced from the original body axes by the trim angle of attack, a1, so that in the new body-axis system the trim attitude angle, 1 = 0 for level flight. We also note that the perturbation, w, in flight-path elevation angle is given by Gw = a + G, where G is the perturbation in attitude angle from its trim value 0 =0. It can then be shown that perturbation equations from level flight take the form: 1 k v +k - (1) g p vv p va w a = vp + ka + q + kze (2) q = a a + kq + kqq + k e (3) qc qa qq 6 q - a (4) where 2CD) k l(D)(5) vv VP L( D1 k (6) kva = CL k 2g (7) V2 kza g = V) (8) Za Vp CL (9) kqa Iy CLe (9) e Vmg CL k - mgc (10) _ mgc2 Cn(1 q a 2IyyVp (11) 2 CMc kqq 2yyVp CL (12) = m gc CNe (13) kqe= yy (3) -H19

Solution of the Linearized Longitudinal Flight Equations Here we have replaced a - qdt by Gw and have added Eqn. 4 to compute 9w from q and a. Using the following data (for the F86-D, flying at Mach 0.8 and 10,000 feet), represent the abov equations on an analog computer and plot the various longitudinal variables for Mach 0.8 and 10,000 feet. How would one excite only the short or long-period pitching motions? Flight constants for the F86-D, flying at Mach 0.8 and 10,000 feet are: a = 0.017 radians k = -23.8 m = 465 slugs kze = -0.238 I = 29,400 slug ft2 k = -0.983 yy x v = 862 feet/sec k = -2.17 p za pS = 0.505 slug/ft kqa = -0.363 c = 8.09 ft kg- = -54.7 gk = -0.0143 kq = -1.57 gkzv = -0.0027 Solution The electronic differential analyzer circuit for solving the above equations is shown in Figure l. Note that the velocity, vp, actually appears as v/g in the circuit. This is done to eliminate the dimensions of feet, so that all variables in the circuit will be dimensionless, or have dimensions of seconds or seconds. Since we are computing in real time (one second in the problem equals one second in the computer solution), this actually means that all dependent variables on the computer can be considered dimensionless. Because of this we can arbitrarily assign a relationship between volts and radians, say 100 volts = 1 radian. Then 100 volts 1 rad/sec, and if vp/g = 100 volts, v/g = 1 or vp = g ft/sec. Note also that angle of attack perturbation appears as 5a in Figure 1. This is done to better distribute the gain between integrators 4 and 6. Computer recordings of the various longitudinal variables for Mach 0.5 and 10,000 feet are shown in Figure 2 following an initial a. Both the short-period and long-period pitching motions are apparent. To excite only the short period, one should start with an initial q; to excite only the long period, one should start with an initial 9w. Computer solutions for Mach 0.M will be similar, except for a longer phugoid period and a higher frequency short-period motion. Revised F86-D Characteristics for Dynamic Stability and Autopilot Studies, Report NA-50-107A, North American Aviation, Inc., January 15, 1952. -H20

Example Problem No. 107 -gkvv -kv~t I I I - - kza/5 -1p I- -5a -' q q - 9~q~~~~ Ot^v q — 0-AA. g. _,'p i Vp.- g V i. I it 400 t4V 0- ^: q:: 0 -::: 0 E — I q- C -- za/5 I;5c 0- - sAAAAA:::;: pAA::A:: 1 A: 1 I I- L-W eo —40'K i::;:; \:,::'::..::::. f i:::/ j.: —- C. ^ —--- IS: - -: O:: 6;;:':::i:::: 0 _ -: t ^:: +5' z: +q o::!'!:,:f ~iguTime, Seconds Figure 2. Typical Longitudinal Transient and Phugoid Response - Mach 0.5 - Altitude 10,000 - C 2 t0. and Phugoid Response - 0.5 - MAlichde 10.5 - Atu 1 - CG 24

Example Problem No. 108 SOLUTION OF THE LINEARIZED LATERAL FLIGHT EQUATIONS by R. M. Howe Problem Statement Let us consider the linearized lateral flight equations as obtained from perturbations from steady-state symmetrical level flight with trim lift coefficient, CL. The dependent variables are roll-rate p, yaw rate, r, and angle of side-slip,p. The equations of motion take the form: I P = + r + k + p + k r k + k a +k ( r = + krp P + krp P + r + a + k (2) IZ r rp rr rcaa r k rr r() = ky: + kyp + kyr r + p w + kS r (3) w= p cos a + r sin al (4) (w is the flight-path axis bank angle. The various k's are given in the following equations: w k C (5) k b (12) xx L zz L _ mgb2 (6) k = mgb CNa (15) pp 2I V C () kr = C (13) xxp L a zz L pr mgbXVP CL (8) k IZZ Cg k (8) k +sP (15) P&a Ixx CL () k V CL - mgb CL k -g CY sin a (16) I a yp C 1 P xx L Lp (6 Lr _mgb Cg gr _00 1 (17) kp'hr - 2x k = (10) 2V p L IzzCL L k =mgb2 CNky Y5r (18) LCL rp 21 V C (1) where the stability derivatives with bars are referred to body axes. Simulate the system on the analog computer using the data given below. How are the "Dutch-roll," spiral, and uncoupled rolling modes excited? -H22

Flight constants for the F86-D flying at Mach 0.8 and 10,000 feet are: a1= 0.017 radians k 22.8 kr3 I = 8630 slug ft2 k = 0.305 I = 35,000 slug ft2 kr = -0.917 m = 465 slugs ky = -0.593 V1 = 862 ft/sec k = 0.017 pS = 0.505 slug/ft k = -1 syr b = 37.1 ft kp9 = -7.42 kpp = -66.0 kpgr= 1.72 kpp = -7.60 kr6= -0.0702 kpr = 1.84 kr6= -0.95 ky = 0.061 Solution The electronic differential analyzer circuit for solving these equations is shown in Figure 1. All of the computer outputs are angles, angular rates, or angular accelerations, and hence, if the problem is solved in real time (one second in the computer solution equals one second in the problem), it is only necessary to choose a single relationship between computer volts and radians, say 100 volts = 1 radian = 1 radian/sec = 1 radian/sec. Note that the roll-rate p appears as p/5. This is because the roll-rate will in general be much higher than the yaw rate, and at the same time it is convenient to have the outputs of all amplifiers at roughly the same voltage level. For the same reason 5P and not P appears as an output voltage. To excite the lateral oscillatory (Dutch-roll) mode one can start the differential analyzer solution with an initial yaw-rate r. To excite the spiral mode it is necessary to start the solution with an initial bank angle 0w. To excite the uncoupled rolling motion the computer is started with an initial roll-rate p. Computer recordings of the various lateral variables for Mach 0.5 and 10,000 feet are shown in Figure 2 following an initial bank angle Ow. Solutions for Mach 0.8 will be similar except for a higher frequency in the Dutch-roll motion. -H23

Solution of the Linearized Lateral Flight Equations Ixz /20xx -k,,j/100 -5/3 -5^ o —-(> ^ - ^ 0.25 - kp/16 |0.25 -p/5 p P kpr/20, 20 5 r 0 20 I kpa8/20 sa 8r 20Ixz/Izz, -p/20 krp /50 0.1 -53 o —--, kp/2 p01 VW -p/5 I s -krr - 4 5 > Br'5, p/5 5' Solving Sin L r I Flgh E.r o-()-1W >1r -H 59, 1 5, - 5ky~ I 8r Figure 1. Electronic Differential Analyzer Circuit for Solving the Linearized Lateral Flight Equations -H24

Exa-mple Problem No. 108 Simulated Lateral Motion of F86-D Showing Short Period Response and Decaying Spiral Mach No. 0.5 - Altitude 10,000 - CG 24% ll.;:.l,^^:lfa~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~lllZll^^^:, I..,..^.^,.. 20 30 140 0 80 Figure 2. Typical Lateral Transient -H25

Example Problem No. 109 TWO-DIMENSIONAL RE-ENTRY TRAJECTORIES FOR LIFTING VEHICLES by R. M. Howe Problem Statement In this problem we are concerned with the calculation of the trajectory of a lifting vehicle as it re-enters the earth's atmosphere from satellite velocity. The computation of vertical forces acting on a satellite requires calculation of the small difference between centrifugal acceleration and gravity, both of which involve nonlinear terms. The usual analog mechanization may yield unsatisfactory results for this reason. However, by using a modified flight-path axis system with suitable equations one can overcome this difficulty.* In this system the xh axis is horizontal and in the plane of the motion (positive forward) and the Zh axis is pointed toward the center' of the earth. The motion of the vehicle can then be described in terms of horizontal velocity Uh along the xh axis and vertical velocity Wh along the zh axis. If we let Xh and Zh be the external vehicle forces not including central force-field gravity, and r be radial distance from the center of the earth, then the following equations can be derived: 1 h d + rUh]t=O (1) h df + -- (1) m r ggr0 1jh Zh Wh r2 Uh (2) h r m where m is the vehicle mass and go is the central force-field gravity acceleration at a fixed radial distance r0. For computer scaling purposes it is usually best to choose r0 equal to the mean radial distance of the trajectory. Equation 1 is essentially the angular momentum integral; when the external force Xh equals zero, it makes the horizontal velocity Uh an algebraic function of radial distance r and avoids an open-ended integration. Equation 2 includes the gravity, centrifugal, and vertical acceleration due to external forces. The equations are more convenient for computer scaling if we define a dimensionless radial perturbation variable Sp given by: r-rO 6p = (3) 0 Also, let us define a dimensionless horizontal velocity perturbation, Suhgiven by: d u h- = U h o - (4) Fogarty, L. E., and R. M. Howe, Analog Computer Solution of the Orbital Flight Equations, IRE Transactions on Electronic Computers, Vol. EC-11, p. 555, August, 1962. -H26

Here Uh0 is the velocity of a satellite in a circular orbit at radial distance rO. Finally, we introduce a dimensionless time i' given by IT =-A t (5) r O In terms of these dimensionless variables, Equations 1 and 3 become / 6U t p (l+&p) m d' + (T+06uh (6) h d2~ p mg 0 1+6 dwh 2d 2Su 6u2 Z h h (u) h _ h + h (7) d-r - 2 - (l+Sp)2 1 + p mgo Here wh is dimensionless vertical velocity, i.e., wh = Wh/UhO. At this point our equations are still exact. For purposes of re-entry simulation we will now assume cp< l< (in typical re-entry trajectories Sp ranges over +0.02), whereupon Equations 6 and 7 become: Xh uh X -d' -p +u (8) 0 mg0 dwh 2d2 Zh h = - = - 2uh - (Suh) + (9) d r dr- mgh where &uh0 is a constant which takes on the appropriate value to make Nuh(0) correct. Note that Suh in Equation 8 will be approximately zero in orbit and will equal approximately -1 near the end of re-entry. Let us turn next to the computation of the external forces Xh and Zh. Since the flightpath angle of a lifting re-entry vehicle will remain quite small except near the termination of trajectory, where the range is influenced very little, we will assume that the aerodynamic drag force D acts horizontally and that the lift force L acts vertically. Thus we write Xh D -h~~~~~~~ -=~~~~~~ ~ ~(10) mgo mg Zh = (L (11)h (11) mg D) mg0 If we assume an exponential atmospheric model, then the drag acceleration -Xh/mg0 can be written: r0 h0 ('P+ P* ) Xh POe 2 U 2 CDS mg0 2 Va 0 mg0 where po* is a constant, ho is the scale height of the atmosphere, p0 is the atmospheric density when Sp +pO* = 0, CD is the drag coefficient, S is the characteristic area on which CD is based, and va is the dimensionless total velocity given approximately by va 1 + uh (13) -H27

Two-Dimensional Re-Entry Trajectories for Lifting Vehicles This approximation assumes that the flight-path,,remains small enough that cos 1. Simulate this system on the analog computer and, assuming L/D = 2.0, plot several re-entry trajectories. Let CDA/mgO =.0066. Solution The electronic differential analyzer circuit used to solve the translational equations of motion in the plane of the orbit, i.e., Equations 8 and 9, is shown in Figure 1. The only nonlinear operation, that of squaring ~uh, is accomplished using one-half of a quarter-square multiplier. The circuit is set up to operate with 1 second of computer time equal to one unit of dimensionless time. By switching the integrators onto repetitive operation this will be speeded up by a factor of 100, since the 1 mfd feedback capacitors are replaced by 0.01 mfd capacitors. At this speed a 10,000 mile re-entry trajectory will require roughly 30 milliseconds. The circuit is scaled with the reference voltage (100 volts) considered equal to unity. Thus by computing 806p we allow Sp to range through +1/80, equivalent to +50 statute miles. Pots A2, A5 and A6 set initial values of horizontal velocity, vertical velocity, and altitude, respectively. The most challenging of the equations presented in the previous section for solving the re-entry trajectory is Equation 12 for the aerodynamic drag acceleration, since this involves the multiplication of an exponential function which ranges over a number of decades with a velocity squared term which also ranges over several decades. A reasonable solution is to utilize fixed logarithm function generators. Thus from Equation 12 Xh pUh2 CDS r0 Xlog ( —-) = lo - 0.434 (p ) log v ( ) mg0 2mg00 a0 a If we choose p0* such that po, the density when Sp + pO0 = 0, is given by 2mg0 U CDS (15) %0 then Equation 14 becomes Xh r0 log(- - 0) = ~434 - ( 6p+ po*) + 2 log va (16) mg0 u0 a where po* is proportional to log (CDS/mgo). The electronic differential analyzer circuit is shown in Figure 2 and employs fixed diode logarithm generators accurate to better than 0.2 percent over two input decades. Amplifier A7 computes total dimensionless velocity, va, in accordance with Equation 13. The two diodes prevent the amplifier output from ever going negative. Pot A7 sets 80p0* and hence represents a quantity proportional to -H28

Example Problem No. 109 0. I I. [ I Xh I 0. 0.1.x2 +10 ~,(6u.)2. m~:h r: - - ~ 90 T-dh -- - Sh+25wh 80log25 0.mg) For the simulation illustrated in this problem we will assume that the drag 0.1Suh0 Xh 0.5 L/D. 1 +100 2 9 rngo 0.0625 1 0.1 0.1 0.1 0.8 aerodynamic forces outside the atmosphere. 0x= u (1 6. (17) 4oh, i-40(&P_6e,))80 0. +100 + 100.uh9. dT -$ + SUh dTh -S -2Suh -(Suh) +m Figure 1. Circuit for Solving the Translational Equations in the Orbital Plane log (CDS/mg0). For the simulation illustrated in this problem we will assume that the drag coefficient CD remains constant, representing an approximately constant re-entry angle of attack. Note that a Mach-dependent CD would require only an additive function of Mach number into amplifier A4. Pot A8 provides a convenient method of obtaining the required gain of 382/80 in summing 80(Sp + p0*) into amplifier B5. By using the log generator in the amplifier feedback loop the output becomes -Xh/mg0 as required. The diode prevents the amplifier output from going more than about one-half volt negative, whereas the two diodes used with amplifier B6 prevent Xh/mg0, from ever going positive (i.e., the drag -Xh/mg0 from ever going negative). When this tends to happen (as it will when the range of the log generator in the feedback of amplifier B5 is exceeded for altitudes that are too high) the lower diode around amplifier B6 conducts and the point labeled Xh/mgO is simply connected to the amplifier summing junction through a 100K resistor. Thus Xh/mgo becomes less than 100 microvolts, which is important to insure negligible aerodynamic forces outside the atmosphere. The dimensionless downrange horizontal distance traveled, s x, is given approximately by SX = j0 (l + ) d(17) -H29

Two-Dimensional Re-Entry Trajectories for Lifting Vehicles 0.1 0.1 0.1 0.1 rogX 0.1'~-100x 1 01 0.1 v 0.1 Xh X -8 \Auh a -0log lOv v3- ( h+ 21.OK 8e ~0. 1 0.1 (Ungrounded Pot) -808P., - - +80P* V 1 + 6uh 1100 0.1 80((8?+eo,*) 10 X +100- 80(S~~~olog -. —--- = 2 log 10va-382 (6$ + (o*) \ mgo / I~log gCDS) mg Figure 2. Circuit for Computing Drag Acceleration Consider first the solution for the zero-lift, zero-drag case, i.e., Xh=Zh=0. Further, assume that the vehicle is in a near-circular satellite orbit, in which case Suh <1l. Neglecting (guh)2 in Equation 9 and eliminating 9uh between Equations 8 and 9, we obtain 2- + sp = 2uh= (18) Thus, for small Suh and Xh=0, the circuit of Figure 1 should behave like a simple harmonic oscillator with natural frequency of one radian per second. The resulting 6p represents the periodic deviation of satellite altitude from circular reference orbit as the satellite goes through apogee to perigee to apogee. Next, consider solutions for finite lift and drag. A recording showing several typical re-entry trajectories for L/D = 2.0 is shown in Figure 3. It is also interesting to record and display other variables, such as vertical velocity, wh, velocity, va, etc. For the case of a non-lifting re-entry vehicle, the flight-path angle g does not remain small. In this case Equations 10 and 11 become Xh D OS D (1~uh) h = - D cos ~ = - mgO -- (19) mg0 mg0 mg0 va and Zh D=- D sin X D h (20) mgO mg0 mgo va Xh Equation 16 for log (-m ) can then be written log (- m-g) = 0.434 h0 (p + PO*) + log va + log (1 + 6uh) (21) -H30

Example Problem No. 109 In the same way we can write Z r (22) log (- mg) = 0.434 (p + po*) + log va + log wh (22) Thus by means of three additional fixed log generators to compute log (l+4uh), log wh, and log-1 mh ), circuits similar to Figure 2 can be used to compute re-entry trajectories for non-lifting vehicles.* Total velocity, v a,must be computed from the equation Va = (l+uh)2 + w2 (23) which can be implemented by fixed diode squarers. Since the peak re-entry accelerations are appreciably higher in this case, appropriate rescaling is necessary. 80( H I ^ "i Sw bilai W X Reentry Trajectories from 338,00tt t 0 Feet T L l =l 2.0, C o = 0.0 066S1 |itrt 1 ii i | 64 2 X t 1 M ~ t S l i | rr 4 t, 10 12 14*4 I 16% 18 20 212 24- tt4_l 72 Iff M s4 Fogarty, L;|-|. E. and R. M. Howe, "Ana log Simulation of the Re-Ent ry of a Ballistic Missile 1 -H31-I lfil L 1 i l Mj 0 - | i~ti iji-vii:ffi -Tigi'lnXW le i i t I i it; i j i* [i iliit! iit", 4 tq. 1 I t ti I -i ~0 2 4 6 - 1- - ~8 100i-lt 2 14 6t182tt0 2 24ttilHti; "t 56;tl ti litiiiFi 31: 3.. Fogarty, Lt E and R MHowe,'Analog Simulation of the Re-EntryofaBallisrttic Mitfltt e 8k-'-, igstldjtll~illil ll~lltl;! jtH i l: it 1 -1 I Warhead~ andI Mfult i pl iDeoyis, Proeedin o the SprEin Int Copue Cf ie, It Ie-, 19 6 40 w ~ ~ ~ 0 I MT-_ 0t: 1 i1 "Ii Reentry Trjectoriesfrom 338,00 Feet LD~e.O, Co/m++J 4+4f 4 1,11f 3 2~~~~~~~~~-H1

Example Problem No. 110 FREE VIBRATION CHARACTERISTICS OF A UNIFORM CANTILEVER BEAM WITH ELASTIC ROOT RESTRAINT by G. Isakson Course: Structural Dynamics Credit Hours: 3 Level: Senior Elective - Graduate Problem Statement Using Targoff's matrix formulation (Jour. Aero. Sciences, October 1947) of the HolzerMyklestad method, set up a program to determine the natural frequencies and mode shapes of vibration of a uniform cantilever beam with elastic bending restraint at the root. Divide the beam into a variable number of equal segments and distribute the mass of each segment equally to its ends. Determine the variation of the frequencies and mode shapes with the magnitude of the spring constant of the root restraint. Express all parameters and variables in nondimensional form. Solution The uniform beam,as shown in Figure la, is divided into n equal segments and the mass of each segment distributed equally to its ends to yield the lumped mass representation shown in Figure lb. L ------— 1 Figure la. Uniform Beam - I\ -i -0 —0i+l i Figure lb. Lumped Mass Model Figure 1. Beam Representation The lumped masses are thus all equal, except for that at the tip, which is half as large as the others. The following linear relationship between variables at stations i and i+l may be written for the general case of a non-uniform beam, -H32

A =i+l [Ei F]i ai (1) where 2/EIol {Aj = ML/EIo (2) I/L 1 o o 0A2 0 1 0 0 [F]i= (3) 0 0 1 0 -0 0 0 1 1 0 0 0 [E], = (4) [EL 2A -E -2/EL 1 0 /ffi EI f/2EI -it 1 and V = shear M = bending moment'1 = slope S = deflection L = length of beam X = length of beam segment X = A/L p = mass of beam per unit length at station i P0 = mass of beam per unit length at root P = P/PO EI = bending stiffness at station i El0 = bending stiffness at root EI = EI/EIo ca = natural frequency, radians/second A: 0X ~/pOL El0 Combining matrices [E] and [F], Equation 1 may be written as follows {t1i+l [i {^1 (5) where 1 0 0 -A2 1 0 2 /2EI -J/EI 1 -_ X /EI 63/ 6E I /2E -J 1-+j p 6EI Applying the condition V=M=O just outboard of the tip mass, we can write the relation, {lar ( i-2 Ali) (7) -H33

Vibration Characteristics of a Uniform Cantilever Beam where subscript r denotes root, and subscript t tip, and 0 ipA/2] 2-A2/2 / [B]= ~0 1pA2 (8) L 1 - pA /12EIj Specializing again to the case of a uniform beam, we can set p = EI = 1 in Equations 6 and 8, and Equation 7 becomes {}r -C C t (9) where [C] = [A]"-1 [B] and the subscript has been dropped from the [A] matrix. Specializing further to the case of elastic bending restraint at the root, with spring constant K, we have Mr = K6' (10) and, from Equation 9, eliminating the first row and setting r=0, iKL/EI C2l C22 r t/'r 1 = C31 C32 t L (11) 0 J C41 C42 Thus, = - 31 C32 r ~1 [oC C2 ] tZt /LA r or 61 t/S r "1 C 1 1 {c4 1 p042 C31C42-C32C41 C31J and KL/EIo = [C21 C22 t/Lr C42C2 -C41 2 (12) L =.C2 22 C42C31-~41C32 In the special case of rigid root restraint, K = oo, and we have the condition, C42C31-C41C32=0 (13) Similarly, in the special case of a hinged root with no bending restraint, K=O, and we have the condition, C4021 -4C1022=0 (14) -H34

Example Problem No. 110 The left-hand side of Equations 13 and 14 may be termed the "residual." In these special cases, the determination of the natural frequencies involves a search for the zeros of the residual by a trial procedure. Following determination of the frequencies, the mode shapes, normalized so that &t/L=l, may be determined by extracting the following relation from Equation 11, C4't + C42t/L 0 (15) Thus C42 (16) t - (C41 I{At is now known and fAji for i=l,2,...,n may be determined by applying successive matrix multiplications as indicated in Equation 5, yielding the non-dimensional values of shear, bending moment, slope and deflection at the selected stations along the beam. The method is implemented on the IBM 704 computer by means of a MAD program. Two main programs are used, one of which determines the value of the spring constant of the root restraint and the residuals for fixed and hinged root conditions for trial values of the frequency parameter, and the other of which determines the mode shape. Both of them make use of two external functions, SETM1. and SETM2., which set elements into the matrices [A] and [B) and which are defined for the general case of a non-uniform beam, so that parameters p and EI are used. The first of the main programs also calls upon an external function RESUN., which performs the cumulative matrix multiplication. The main program for spring constant and residual determination will accept either an assortment of trial values of the frequency parameter or an initial trial value to be incremented by a specified amount a given number of times. The external function SETMI. sets those elements of [A] and [B] which are independent of the frequency parameter A. The external function SETM2. sets the remaining elements of A] and IB]. The following additional symbols used in the flow charts and programs are defined: CRITER = Switch determining whether frequency data list is complete or is to be incremented. M = Number of frequency values to be read in. RCANT = Residual for cantilever beam with rigid root restraint. RHINGE = Residual for beam with freely hinged root. SPRING = Spring constant of root restraint in non-dimensional form. NN = Number of mode shapes to be determined. CASE = Identifies by means of the alphabetic names CANTIL, HINGED, SPRING, whether the root of the beam is rigidly restrained, freely hinged or elastically restrained respectively. -H35

Vibration Characteristics of a Uniform Cantilever Beam SHAPE(I,1) = VL2/EI at station i. SHAPE(I,2) = ML/EIO at station i. SHAPE(I,3) = g' at station i. SHAPE(I,4) = G at station i. Flow Diagram - Main Program for Frequency Determination |READ | PR5gr T RVIr - T C Ni MI N? \A m _ y EX/rCcure /3ZP It EX6EC UTE <:< SCT fA.1, C rM2 ExECuri C -- C ( ps I. r AINT Es QLTS 2 SPpptN CcKT'T4Nf A)r cPSir C C oNsr6/r5 -N' Sr1 iPE, lDGUALS A /i I 4L5, C4 2 CL - 2 Flow Diagram - Main Program for Mode Shape Determination RCA _ o I j- ExECUTI /ZCr1\ C P FmrfT "- A PCI,, CA5C, S1/ jrv t, Co 1 ) ^,'N N SE \D / iM C / y L- I N K _ _ N1_ _ P NN C-41 / 42 PPr/fr I Ex c< r7 CPc vTr / rT7TA\ l-^CA 5 F5pO Ccrj4 LO. - C TtM' A PIEO T) I0Z 2.... 2, 5 \1- A,4 T- - + 5I AJ, 2 \L rK \= |i' 3,J. p^Er CCL/V2^\ _ RA __ K / — VI/ r ^ /ib/N / S |. 5H.4 -H36

Example Problem No. 110 Flow Diagram - Subroutine RESUN. fT/A F/ B NT( j- I7 1) CT j77v sr4 5Z2AMA/I\A DEL-TA /ELT\ C1~y 1 Y 111 1) 1 $O. &2: LPSIL. 5 LS 2(3KK 4 1) I~K A11K CK EIT L D Sr MLcELTA MAD Program - Main Program for Frequency Determination $ COMP.ILE MAD, EXECUTEt DUMP, PRINT OBJECT RPROGRAM FOR UNIFORM BEAM NATURAL FREQUENCIES BY RHOLZER-MYKLESTAD METHOD DIMENSION A(16,DESC),B( 8DESC(3)),C(8,OESC(.3) )LAi'MDA(200) VECTOR VALUES DESC=2,1,4,2,1,2 INTEGER N,, MQ,CRITER START READ FORMAT INPUT1,N,CRITER,M VECTOR VALUES INPUT1=$I2,I1 I3*$ PRINT FORMAT OUTPT1,N,M VECTOR VALUES OUTPT1=$63H1NATURAL FREQUENCIES OF UNIFORM BEAM 1 BY HOLZER-MYKLESTAD METHOD/4H N =,I3,S53,3HM =,I4*-$ WHENEVER CRITERE,.1 READ FORMAT INPUT2,LAMDA(1).**LAMvDA(M) VECTOR VALUES INPUT2= (9F8.5)*$ OTHERWISE READ FORMAT INPUT3,LAMDA(1) DELLAM VECTOR VALUES INPUT3=$2F6.5*$ THROUGH PSI, FOR Q=2,1,Q.G.M PSI LAMDA(Q)=LAMDA(Q-1)+DLLLAM END OF CONDITIONAL LBAR=1 /N EXECUTE SETM1i(LBAR,1Z,1,o, AB) THROUGH ALPHA, FOR Q=1,1,Q.G.M LAMDAA = LAMDA(Q) EXECUTE SETM2.(LBAR,1,1,,LAMiDAAA B) EXECUTE RESUN.(A,B,C,N) RCANT=C(4,2) C (3,1)-C(4,1 )*C(322) RHINGE=C(4,2 ) C ( 2 1)-C (41 )C( 2 2) SPR I NG=RH I NGt/RCANT ALPHA PRINT FORMAT OUTPT2,LAMDAA,IRCANT,ItHIiGGa,SPi-<ING,C(4,1),C(4,2) VECTOR VALUES OUJTPT2=SSHuLAMDA =,F9,3/11H REsI DUALS-',j, 112HCANTILEVER =,E1,, iS5,HHINGED =,tE15.lb/ ii-i.SPI<INo CFNST' Ai',T 2 =,E15.8/9H C(4,1) =,E15.d,S1u,6HC(4,2) =z,El1.8'$ TRANSFER TO START END OF PROGRAM -H37

Vibration Characteristics of a Uniform Cantilever Beam MAD Program - Main Program for Mode Shape Determination $ COMPILE MAD, EXECUTE, DUMP, PRINT OBJECT RPROGRAM FOR DETERMINING UNIFORM-BEAM MOD SHAPES BY RHOLZER-MYKLESTAD METHOD DIMENSION A(16,DESC),B(8,DESC(3)),SHAPE(85,OESC(6)) VECTOR VALUES DESC=2,142,2,1,2 24,4 INTEGER PI,J,K,N,NN START READ FORMAT INPUT1,N,NN VECTOR VALUES INPUT1=$I2,I3*$ LBAR=1./N EXECUTE SETM1.(LBAR,1.,1, 0,A,AB) THROUGH ZETA, FOR P=1,1,P.G.NN READ FORMAT INPUT2,CASE,SPRING,LAMDAA,C41,C42 VECTOR VALUES INPUT2=$C6,E15.8,F85,2E15*8*$ PRINT FORMAT OUTPT1,CASESPRING,N,LAMDAAC41 C42 VECTOR VALUES OUTPT1=$1HO, C6 S5,1 7HSPRI NG CONSTANT =,E16.,S5 1,3HN =,I3/8H LAMDA =,F10.5,S5,5HC41 =,E16.8,S55,HC42 =,E16.d2$ EXECUTE SETM2.(LBAR,.,1.,*LAMDAAAB) SHAPE (01) =0. SHAPE(0,2)=0. SHAPE(03) =-C42/C41 SHAPE (0,4)=1 SHAPE( 1,1)=B( 12) SHAPE(1,2) =B(2 2) SHAPE ( 1 3) =SHAPE ( 3 ) +B(3,2) SHAPE(1,4)=B(4,1)*SHAPE(U,3)+B(4,2) THROUGH THETA, FOR J=1,1,J.G.N-1 THROUGH THETA, FOR I=1,l,I.G.4 S=O. THROUGH TAU, FOR K=1,1iKG.4 TAU S=S+A(IK)*SHAPE(JK) THETA 5HAPE(J+1I )=S PRINT FORMAT HEAD VECTOR VALUES HEAD=$1HSu5,7HSTATION,S,8 5HSHEAR, 51 l, 6Hi-;O.;Ei'T, 1S9,5HSLOPE,S10 1HDEFLECTION*$ THROUGH ZETA, FOR I=u,1,I.G.N ZETA PRINT FORMAT DATA,I,SHAPE(I,1 )*..SHAPE(it4) VECTOR VALUES DATA=$1H, 9 F18 5,2 F15 5, F 1 7 * TRANSFER TO START END OF PROGRAM MAD Program - Subroutine RESUN. $COMPILE MAD, PRINT OBJECT RSUBROUTINE FOR UNIFORM BEAM FREQUEiNCIES Y H-iLZER-iYKLESTAD RMETHOD. MATRIX PRODUCT FOR RESIDUAL DETERii,'iIilATION. EXTERNAL FUNCTION (AB,CN) DIMENSION D(8,DES) VECTOR VALUES DES=2,192 ENTRY TO RESUNs THROUGH BETA, FOR I=1=1,IG.4 THROUGH BETA, FOR J=1,1,J.G.2 BETA C(I,J)=B(I J) THROUGH GAMMA, FOR P=1,1,P.G.N-1 THROUGH DELTA, FOR I=1,1,I.G.4 THROUGH DELTA, FOR J=,1,J.G.2 S=O0 THROUGH EPSIL, FOR K=1,1,K.G.4 EPSIL S=S+A( I,K)'C(K,J) DELTA D(I,J)=S THROUGH PHI, FOR I=1,1I*.G.4 THROUGH PHI, FOR J=1,1,J.G.2 PHI C(I,J)=D(IJ) GAMMA CONTINUE FUNCTION RETURN INTEGER I,JP,K,N END OF FUNCTION -H38

Example Problem No. 110 MAD Program - Subroutine SETM1. and SETM2. $COMPILE MAD, PRINT OBJECT RSUBROUTINE FOR BEAM NATURAL VIBRATIONS BY HOLZER-MYKLESTAD RMETHOD. SETTING ELEMENTS INTO A AND B MATRICES EXTERNAL FUNCTION (LBAR EI B ARROBAR tLAML)AAA B) ENTRY TO SETM1. A(1I1)=1. A( 2, 1 )=LBAR A( 12)=0. A(2 92)=1l A( 32) =-LBAR/EIBAR A( 3,1 )=.5*LBAR*A( 3.2) A(42 )=-A( 3,1) A(4 1)=(1*./3 )*LBAR*A( 42) A( 1 3)=O. A(2.3)=0. A(3,3)=1l A(43 )=-LBAR B( 1 1)=0O B(2,1)=0O B(3,1)=1. B(41 )=-LBAR FUNCTION RETURN ENTRY TO SETM2. LAMS=LAMDAA P 2 A( 1,4 )=LBAR*ROBAR*LAMS A(2,4) LBAR*A (1 4) A (34)=-5*( LBAR/E I BAR)*A (24) A(4,4)=-( 1/3 )*LBAR*A(3,4)+1. B( 1 2) = 5*A( 1,4) B( 22 )= 5*A(2 24) B(3s2)=*5*A(394) B(4,2)=5*A( 44 )+,5 FUNCTION RETURN END OF FUNCTION Computer Output - Typical Results from Mode Shape Determination Program CANT IL SPRING CONSTANT = -041597L 36E 06 N = 10 LAMDA = 3.50000 C41 = -0.;10077o4E 01 C42 = 015194210E 01 STATION S-i EAR M"iOMENT SL OPE DEFLECTION 0.0.00000 u u000 1.38032 1.00000 1. 61250 u-.06125 1.37725 0.86207 2 1 66854. 22810 i. 3278 0.72493 3 2.55658 u.4-6376 1.32719 0.59022 4 3.27959 0. 1172 1 26242 0.46046 5 3.84366 1.19609 1.16203 0.33892 6 4 25884 1.62197 1.02112 0.22941 7 4.53986 2.07596 0. 3623 O.13616 8 4*70666 2.54662 *0.o510 0.0637u 9 4.78470 3.02509 O.J2651 0.01672 10 4 80519 3.50561 -0.U0002 -0 CQOOO -H39

Vibration Characteristics of a Uniform Cantilever Beam Computer Output - Typical Results from Frequency Determination Program i3' r 1RHL F.EUEHC I E F H F. i:OFi i 3. H O iLZ E I;, 7 ET H, = 1. 00! -_ ES iD 1_L', s-.-. I.' I L, E r = _. i':'7' 5 5 6 E'-; 5 G E [ j 3. 3-;.5,:7 i. i'. F ".: G.' 0 H i- T R T. 1.-., 1,, 3 2544 7E 0 7 i_:!;:, ~.....O. 1, 1 0 i"'7 7.1i'i! C" C,::i t. 2:' - 0' i - 1 9 j-0 9 2 E i i. 1 I.:, E? - i T E ii =, 5 i::. H IT C.D....... i 7 S i 1 i Ei 0 i. S F: I.H G 1: _ t T R 1-1 T.-. 4 r 9:,'.,:i E i 3 LS R ].i R - 2 SRT 01.&7 S.0 r / E S D i S - i:'' I L E E.F, E - 1 I E I -..: i 7? 3 S PRIG COHST RT.. 0.41-525.3E. 04 i, 1 - — 0. j 1.-. Ci i i::. 2 -: 2,'; 9.',,_, j E i';'NL!:rhi4 i- 1.... iR E -; n D I R L. _-F I -.' i L E',L,F E...: E R -'.8255 -i! 0 1 E'!' i'! Ei E!: - - 7' i; L' L -;. 7 i E.; 0 0 -. E S- I D U. RA- I f'i L' I: -'E - 1 E 0 H - -. E 0'R! s-... 1 Gr:.' C S T 1.- i,: f'-,'4I::, -..F' JT1:I, T - 0 31.}.::0iE 02 |sj 1, = -0 = 4. i L.; Di' i i 1i 5. i. S - i f L' - 1 I' * E 01 I i E: S T E 0 I I i l f!. 1!J.- L'.I I - r- 1 I' G ". 0 i N S= T -i. E.i -H40 -:z ~- If a " i ~~~~~~~~ I~~i r~~~nt i'-, I::::~~~~~~~ ": - ~: I 1A~ ~ ~~~~~H0

Example Problem No. 110 Discussion of Results The spring constant data are shown plotted versus frequency parameter in Figure 2. It is seen as expected that a beam with positive elastic restraint in bending at the root will have an ith mode frequency which lies between the corresponding frequencies for the hinged and rigidly cantilevered beams. The residual data were used to determine the first two natural frequencies for the hinged and ridigly cantilevered beams. A comparison between these results and exact analytical results is shown below: Frequency Parameter Root Condition Mode Present Exact % Error Method Solution Hinged 1 0 0 0 2 15.219 15.418 -1.32 Rigidly Restrained 1 3.500 3.516 -0.46 2 21.690 22.030 -1.55 Mode shape data were computed only for one case, namely, the rigidly cantilevered beam. The data agree with the exact solution within plotting accuracy.'~ ~ it i~F-~fift ~ m_-f'+KI J~]!!~~'i'IiI''p''I ~" l —l-!' lHl. am'~m-'X~ g t i' I' t g -' XH I ~ - X 4 - < Ii: ~ m t Ad(4 ~ C aO C >O o o H4-H+H-H4 u' - rl.' ~..eF" j [, ^'t i fi rrr-~rl-"rt1~-1"inil.tltM TFl ftli i fim-ti/illi Ad ~ 0 1.0 2.0 3.0 15 16 17 18 19 20 21 Frequency Parameter, Figure 2. Spring Constant Variation with Frequency -H4l

Vibration Characteristics of a Uniform Cantilever Beam Comments and Critique The method treated in the present example is gaining wide acceptance in practice and, in a more general form, is applicable to a wide variety of static and dynamic structural problems. It illustrates to the student how a method which would be very laborious for hand computation is readily implemented on a digital computer and can yield a large volume of data in a very short tii Its extension to the more general case of a non-uniform cantilever beam and beams with othe: support conditions can be readily accomplished. By varying the number of segments into which the beam is divided, the student can determine the effect of this number on the accuracy with which specific mode frequencies and shapes can be determined. The present example could be used appropriately in any course concerned with vibration theory for continuous elastic bodies. -H42

Example Problem No. 111 RE-ENTRY FLIGHT PATH OF A LIFTING VEHICLE by Gabriel Isakson and Joseph L. Lemay Course: Flight Mechanics of Space Vehicles Credit Hours: 3 Level: Senior or First Year Graduate Problem Statement Given the values of altitude, velocity, flight path angle and direction, and geocentric latitude at the initiation of re-entry of a lifting vehicle into the atmosphere, determine the re-entry flight path. Include the geometric effects of earth oblateness, but neglect oblateness gravity effects. Assume the ratio of lift to drag of the vehicle to remain constant. Consider only symmetrical motions and assume that the vehicle is at all times trimmed so that there is zero moment about the pitch axis. Mathematical Formulation of Problem The coordinates and motion variables of the vehicle are shown in Figure 1. — LOCAL HORIZONTAL / M(NORMAL TO RADIUS VECTOR) FLIGHT V PATH A INITIAL RADIUS /r \?R / VECTOR EARTH SURFACE EARTH CENTER Figure 1. Coordinates and Motion Variables of Vehicle The following two kinetic equations for forces normal and tangential to the flight path may be written, L + m V (~ + 2) - mg cos = 0 (1) -D + mg sin - m V = o (2) where L = aerodynamic lift D = aerodynamic drag m = vehicle mass V = vehicle velocity g = acceleration due to gravity -H43

Re-Entry Flight Path of a Lifting Vehicle 1 = WNODAL: \LINE \UATOR * ^ ^/ ASCENDING NODE Figure 2. Geometric Parameters for an Inclined Orbit and dots denote differentiation with respect to time. In addition, the following three relations are needed, r = -V sin & (3) rG = V cos Y (4) r - R + h (5) as well as the aerodynamic relations, L pV2AC = pV2 ACD(D) (6) D = pV2AC (7) where A = reference area for CL and CD CL = lift coefficient CD = drag coefficient p = air density The variation of air density with altitude may be approximated as follows through most of the altitude range of interest, = h e (8) -H44

Example Problem No. 111 where P0 = 0.0027 slug/ft3 and P = 1/23,500 ft1. The radial distance R from the center to the surface of the earth can be related to the geocentric latitude, Y, as follows, R = RO (l-fsin 2) (9) where f = 0.003367 R0 = 20,926,428 feet From Figure 2 it can be seen that XL can be related to the angular coordinate 0 in the orbital plane as follows: sin Y = sina sin / (10) permitting the rewriting of Equation 9 in the form: R = Ro (l-fsin2a - sin2/) Furthermore, 0 = ot+ (11) where 00 is the initial value of f. Combining Equations 5, 9, and 10, we have h = r-R0 {1-f sin2' sin2(00+9)} (12) Equations 1 to 4 inclusive and Equations 6, 7, 8 and 12 now comprise the set to be solved. A reduction in the number of dependent variables can be effected by a transformation in which time is replaced by G as the independent variable. Thus, from Equation 4, 0 = -Zcos (13) and, applying this relation, dV dV V cos = (V2), cos (14) V d - r 2r 0d Y, V cost dr' V os dO r (15) r= dr9 r' Vcos (16) where primes denote differentiation with respect to 9. Substituting Equations 6, 7, 8, 14, 15 and 16 into Equations 1 to 4 inclusive, and defining the drag mass parameter, MD = CD Ap /2m, we have L MD e h V2 + V- cos h (1+2) - g cos 0 = 0 (17) D D e -MD e- h V2 + g sin O - (V2)' cos 0 (18) D 2r r' = -r tan (19) The variable V can be put in a non-dimensional form by dividing it by the square of the circular orbital velocity corresponding to radial distance r from the center of the earth. Thus we define -2 v2_ V2r (20) V - _rg W where g = /r2, with J. the gravitational constant neglecting oblateness effects. -H45

Re-Entry Flight Path of a Lifting Vehicle We can now write (72) -2 r2 2 -2 (21) Substitutingt Equations 20 and 21 into Equations 17 and 18, and rearranging, we have 1 L/D MD e-Br'6 7 2 cos (22) (2) = (2-V2) tan - 2MD (23) cos 2 which, together with Equations 12 and 19, comprise the set of equations to be solved. Round-off error in a digital solution can be greatly reduced by replacing Equation 19 by (Sr)' =-(6r+R0) tan, (24) where Xr = r-Ro, (25) and is small in magnitude as compared with r. Equation 12 can now be replaced by h =r + R f sinin sin2(0 + ) (26) Rewriting Equation 25 in the form r = fr + RO, (27) we have finally the set of equations 22, 23, 24, 26 and 27 to be solved. Initial conditions may be selected on the basis of an approximate solution for a non-oscillatory trajectory around a spherical earth and are then given by the expressions, -2 1 (8) 1+- MD e r ID DD 1 1 F -- (1 + MD e-hr) + R f sin a sin 2(00+)) (29) — Y\ + Pr L D MD 0 LMD2 h r for a given initial value of h. If it is desired to examine the characteristics of an oscillatory trajectory, the initial conditions may be varied from the values given above. Solution Equations 22, 23, 24, 26 and 27 of the preceding section are solved on an IBM 709 computer subject to initial conditions given by Equations 28 and 29. The program is written in the MAD language. After introduction of pertinent parameters, including initial values of h and 0, initial values of V2 and X are computed from Equations 28 and 29 respectively. The differential equations 22, 23 and 24, supplemented by Equations 26 and 27, are solved using a fourth-order RungeKutta procedure (Hildebrand, Introduction to Numerical Analysis, McGraw-Hill, p. 237, Eq.(6.16.8)) which is written into the main program. -H46

Example Problem No. 111 The step size in the integration procedure is decreased at specified altitudes to allow for the fact that the wave length of any oscillation which may be present in the trajectory decreases with decrease in altitude. The computation may be terminated by the satisfaction of any one of three conditions. These are a.) Y greater than 89~, b.) V negative, c.) h negative. Values of the variables,, -2V V, r and h are stored in a matrix until data for one hundred successive steps are accumulated, and are then printed out. The following symbols used in the program are defined: ALPHA = a LIMGAM = limiting value of' for terminating condition F = f DELR = LR H= h R= r PHI = 0 V2 = V INC = increment in GAM = B = P V =V LD = L/D THETA = 9 (in radians) MD = MD THETAD = 0 (in degrees) RE = RO Flow Diagram SAk~ ~Dr, j L' /b /23 Prt- t 14) 5f, /i VD 00- H N,' =,,'.A/,H /y)^ - P- < 1 - e ^ ^ ^l J( C/tY.^ <^-^ CJ - ~ /~'~~.I ^,.~ /~...........t \H< J)g ~ ~ ~ ~ YI I I.'.... -' - l |Bj LDPRoc f /C - I 5-'+3~)-' v' sr M \ J -Tc OC 0072-U4.ZC 0$o34 3 53- o3)3 2r ---------- _- Sr.Kt4n j STCN) 7HU7AJ) jUtifcE-KUT(4A 57 T(Nl)z C^4A/v $ ____ VTP:)T PROCE rrE4'5 T I N —15)) srT_ C s) ^ f Additional symbols in the program are used in connection with the Runge-Kutta procedure. -H47

Re-Entry Flight Path of a Lifting Vehicle MAD Program START READ FORMAT IMNUT.ALPHAF,HPHIINC B=1./23500. LD=3. MD=.0003 RE=20926428. LIMGAM=89.9/(180.)*3.14159265 LDMD=LD*MD X=RE*F*SIN.(ALPHA).P.2 DELR=H-X*SIN.(PHI).P.2 R=DELR+RE Y=LDMD*EXP.(-B*H)*R V2=l./(l.+Y) GAM=((l1+Y)*2.)/(LD*(l1/Y+B*R))+(y*B*X*SIN.(2.*PHI))/(1.+Y*R* 1B) PRINT FORMAT CHECK,ALPHAPHI,F,B,LDMDREGAMV2,DELRRH N=1 THROUGH ENDt FOR THETA=09INCTHETA.G.7* WHENEVER H.L. 220000.,INC=.00872664 WHENEVER H.L 170000.*INC=.00436332 WHENEVER GAM.G.LIMGAM.OR.H.L.O. PRINT FORMAT OUTPUT#ST(1).**ST(N-1) TRANSFER TO START END OF CONDITIONAL GAMDN=1*/V2-1.-Y/COS.(GAM) V2DN=(2.-V2)*TAN.(GAM)-(2**Y*V2)/(LD*COS.(GAM)) DLRDN=-(DELR+RE)*TAN.(GAM) GAM51=GAM+GAMDN*INC*.5 V251=V2+V2DN*INC*.5 DELR51=DELR+DLRDN*INC* 5 R51=DELR51+RE H51=DELR51+X*SIN.(THETA+PHI+INC*.5).P.2 Y51=LDMD*EXP.(-B*H51)*R51 GAMD51=l1/V251-1.-Y51/COS.(GAM51) V2D51=(2.-V251)*TAN.(GAM51)-(2**Y51*V251)/(LD*COS.(GAM51)) DLRD51=-(DELR51+RE)*TAN.(GAM51) GAM52=GAM+GAMD51*INC*.5 V252=V2+V2D51*INC*.5 DELR52=DELR+DLRD51*INC*.5 R52=DELR52+RE H52=DELR52+X*SIN.(THETA+PHI+INC*.5).P.2 Y52=LDMD*EXP&(-B*H52)*R52 GAMD52=1./V252-1.-Y52/COS.(GAM52) V2D52=(2.-V252)*TAN.(GAM52)-(2.*Y52*V252)/(LD*COS*(GAM52)) DLRD52=-(DELR52+RE)*TAN.(GAM52) GAM11=GAM+GAMD52*INC V211=V2+V2D52*INC DELR11=DELR+DLRD52*INC R11=DELR11+RE Hll=DELR11+X*SIN.(THETA+PHI+INC).P.2 Y11=LDMD*EXP.(-B*H11)*R11 GAMD11=1./V211-1.-Y11/COS.(GAM11) V2Dll=(2.-V211)*TAN.(GAM11)-(2.*Yll*V211)/(LD*COS.(GAM11)) DLRD11=-(DELR11+RE)*TAN.(GAM1i) GAMD=(GAMDN+2.*GAMD51+2.*GAMD52+GAMD11)/6. V2D=(V2DN+2.*V2D51+2.*V2D52+V2D11)/6. DLRD=(DLRDN+2.*DLRD51+2.*DLRD52+DLRD11)/6. GAM=GAM+GAMD*INC V2=V2+V2D*INC WHENEVER V2.L.O. PRINT FORMAT OUTPUTST(1)...ST(N-7) TRANSFER TO START END OF CONDITIONAL DELR=DELR+DLRD*INC R=DELR+RE H=DELR+X*SIN.(THETA+PHI+INC).P.2 Y=LDMD*EXP.(-B*H)*R V=SQRT. (V2) -H48

Example Problem No. 110 MAD Program, Continued THETAD=(THETA+INC)*57,295780 ST(N)=THETAD ST(N+1)=GAM ST(N+2)=V2 ST(N+3)=V ST(N+4)=R ST(N+5)=H N=N+6 WHENEVER N.G600 PRINT FORMAT OUTPUTJST(1)***ST(N-1) N=1 END END OF CONDITIONAL INTEGER N DIMENSION ST(1000) INTERNAL FUNCTION TAN.(Z)=SIN.(Z)/dOS.(Z) VECTOR VALUES INPUT=$5F14.6*$ VECTOR VALUES OUTPUT=$(6F20.6)*$ VECTOR VALUES CHECK-$(6F20.8)*$ END OF PROGRAM Discussion of Results Three typical trajectories computed by means of the preceding program are shown in Figure 3. All are for a once-around polar orbit (a = 90~). One involves the assumption of a spherical earth (f=0). The remaining two are for an oblate earth, one starting at the equator and the other at a pole. Otherwise, initial conditions are the same in all three cases. The results illustrate the effect of earth oblateness on the trajectory. Oblate, o0 = 0~ 400 F. / /-Spherical w 3,50_ UL. Z 250 <Z Oblate, o= g900 C:) o 200 15o for 150 -- -- -- --- --- -- --- --- -Around -- — Spher cal I-I?~ Oblate, o =O 90~ —"' 50 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 9-DEGREES Figure 3. Comparison of Spherical and Oblate Earth Trajectories for a Once-Around Polar Orbit -H49

Re-Entry Flight Path of a Lifting Vehicle Critique The present example serves to illustrate the power of the digital computer in solving problems which are very difficult to solve by analytical means. In the present case the problem involves the solution of a set of simultaneous nonlinear differential equations. The program demonstrates the implementation of a typical numerical analysis technique on the computer, in this case a particular form of the Runge-Kutta method of solving differential equations. If desired, the program can be simplified by making use of a library subroutine for solution of differential equations. -H50