[37] As in most scientific and engineering work, analytical studies are required in conjunction with experiments to understand the results of aeronautical research and to predict the characteristics of new aircraft or aeronautical systems. Much of the progress in improving the performance and safety of airplanes relied heavily on analytical work, which in turn depended on the availability and understanding of analytical techniques. This section outlines the status of analytical techniques that I had encountered in my studies and that were generally available to aeronautical researchers at the time I started employment with the NACA in 1940.
Studies of the stability and control of airplanes, my field of specialization, relied almost completely on the application of Sir Isaac Newton's laws of dynamics, which in turn led to the need to solve differential equations. In my college courses in mathematics and physics, Newton's laws of motion were among the first subjects studied, and the various means to apply them in different scientific disciplines occupied most of the subsequent curriculum. Mathematics courses, particularly differential and integral calculus, were presented to students as a general preparation for all the courses. These methods were originally developed largely by Newton to solve his own problems. The present notation for differential and integral calculus was originated by Jacob Liebnitz, a contemporary of Newton. A branch of calculus called differential equations occupied a whole term at MIT. Many of the methods that were taught to solve differential equations, however, applied to special forms of equations that did not occur in connection with airplane stability and control. This subject requires the solution of simultaneous linear differential equations with constant coefficients. The standard mathematics curriculum did not give this subject any special emphasis. As a result, it was left for the professors in the aeronautical courses to emphasize the importance of this particular equation. These professors, in my courses at MIT, were practical engineers without a strong mathematical background. As a result, graduates were provided with only a minimal introduction to the branch of mathematics most useful for their subsequent work. I feel that this situation existed in most of the colleges in the United States at that time, with the result that engineers who wanted to go into this field had to spend a great deal of their time in reviewing work done largely by electrical engineers or by some aeronautical researchers in Europe. Many chose to put their emphasis on other areas with which they were more familiar. The result was undoubtedly a slowing down of progress in this field.
The reason that simultaneous differential equations with constant coefficients arose in [38] aeronautical stability work may be explained in somewhat more detail as follows. For each variable that describes motion (the degrees of freedom), a differential equation is set up. In the case of longitudinal motion of a rigid airplane, for example, the variables would be vertical displacement, horizontal displacement, and pitch angle. Newton's third law states that force equals mass times acceleration. Thus for each variable, the aerodynamic or gravity forces due to the motion, which may produce forces or moments proportional to the displacement or velocity, are equated to the mass times the acceleration of that variable. Acceleration is represented by the second derivative of the variable with respect to time, velocity the first derivative, and displacement does not involve a derivative. As a result, a differential equation of the second order is written for each degree of freedom, which results in three simultaneous differential equations. These equations, in the usual formulation, are linear for the following reasons.
The airplane is first considered as flying in a trimmed condition. After a disturbance, the forces on the airplane change as a result of the effects of changes in angle of attack or airspeed on the various components of the airplane. These changes are necessarily small in the normal unstalled range of flight. The stall angle, at least on airplanes of WWII vintage or earlier, was usually about 15 degrees, and most disturbances producing forces within the structural capability of the airplane would be much less than this value. Ideal fluid theory shows that within this range, forces vary nearly linearly with angle of attack. The only factor that would change this condition would be the effect of the boundary layer, but the boundary layer on full-scale airplanes designed for efficient flight is so thin that it causes little effect on the forces. Likewise any change in airspeed caused by a disturbance is likely to be small compared to the initial large value of airspeed. Therefore, though the forces vary as the square of the airspeed, the variation in force for a small change in airspeed may be considered to vary almost linearly with the change in airspeed. By comparison with many mechanical systems, the airplane as a whole does not have any coulomb friction tending to hold it in its trimmed condition. Coulomb friction is the force caused by rubbing two solid surfaces over each other, and from an analytical standpoint, is usually considered to be a constant force that is independent of velocity and opposing the motion. The force caused by aerodynamic effects is in marked contrast to that appearing in mechanical systems. This difference results in the validity of linear equations to describe the motion of an airplane, whereas linear equations are usually a poor representation of mechanical systems.
Solution of the system of linear differential equations with constant coefficients can be accomplished by the classical method given in most books on differential equations. These solutions go back to the work of mathematicians in the earliest days of mathematics and were perhaps first summarized in a paper by Leonhard Euler in 1739. This solution is considered elementary from a mathematical standpoint, inasmuch as it involves elementary functions such as exponentials and sine and cosine functions. From the standpoint of the practical engineer using the equipment available when I came to work at Langley, which consisted of slide rules or mechanical calculators, the solution is very time-consuming. Determining the characteristics of the various modes of motion, such as the phugoid oscillation or short-period mode, requires solving for the roots of a fourth-degree algebraic equation. If the effects of a simple autopilot are included, a sixth-degree equation results. This problem can be solved only by trial and error or by methods of successive approximations. Further lengthy calculations are required to determine the constants giving the amplitudes of each mode for known initial conditions.
Aeronautical engineers confronted with these problems made efforts to devise simpler methods for solving for the roots of [39] equations, as will be described in the subsequent material. At the same time, and usually without the knowledge of aeronautical engineers, electrical engineers had developed different approaches that allowed practical solutions of even more complex systems of linear differential equations. Such solutions were required in the design of electrical circuits, vacuum-tube feedback amplifiers, and transmission lines. One method, first summarized in a book by Hendrik ~ Bode of the Bell Telephone Laboratories, entitled Network Analysis and Feedback Amplifier Design (ref 5.1) exploited the use of sinusoidal forcing functions of various frequencies to the dynamic systems under consideration, from which the stability of the systems could be determined This method was slow to be discovered by the aeronautical engineering profession because of the unfamiliar notation and applications of the electrical engineers. Later, the method was called the frequency-response method and was widely used. A second approach with the general title operational methods was also introduced by the electrical engineers. A British electrical engineer and mathematician named Oliver Heaviside devised a system of operational calculus about 1887. The Heaviside operational calculus was publicized in this country in a book by Vannevar Bush, entitled Operational Circuit Analysis (ref. 5.2) and was introduced to the aeronautical profession by Robert T. Jones in NACA report No. 560, A Simplified Application of the Method of Operators to the Calculation of Disturbed Motions of an Airplane (ref. 5.3). Later, the book by Murray F: Gardner and John L. Barnes, Transients in Linear Systems Studied by the Laplace Transformation (ref. 5.4), described the operational method based on Laplace transforms, which became generally accepted as easier to understand than the Heaviside method. The advantage of operational methods is that solutions for frequently encountered equations and for special inputs such as steps and ramps can be obtained from tabulated precalculated formulas.
A third approach for solving these systems of equations was the use of simulators, generally referred to at the time of their development as differential analyzers. Vannevar Bush developed a mechanical differential analyzer at MIT that I saw about 1934 while I was still in high school. This machine took the equations in an integrated form so that the various terms required integration rather than differentiation. The integration was done by devices known as rolling-wheel integrators. The low torque output of these devices was amplified by winch-type electromechanical servos to drive a large array of shafting and gearing that allowed setting in the correct constants from the equations. The whole equipment required a room about 25 by 60 feet. This machine had excellent accuracy and was used by an MIT student to get some of the first solutions for the motion of airplanes with automatic controls. Later, such simulators were made with servodriven potentiometers to enter the constants and electronic operational amplifiers to perform the integrations. These machines had tremendously increased capabilities, but they in turn became obsolete with the development of electronic digital computers.
Although the airplane as a whole can be described accurately by linear differential equations, many of the subsystems, such as control surfaces or autopilots, involve nonlinear components. For example, control surfaces have static friction, and electronic devices often have switches that give a discontinuous output. In general, nonlinear systems include devices in which the output varies in a nonlinear but continuous manner, such as a crank, or in which the output varies discontinuously, such as a bearing with coulomb friction or an on-off switch. Much of the mathematical analysis of nonlinear systems when I was in college had been concerned with the continuous type of nonlinear systems, though solutions were usually available only for systems capable of representation by special types of differential equations. The only method of analysis for discontinuous systems with which I was [40] familiar was the phase plane method, a graphical method described in the book by Nicolai Minorsky, Introduction to Non Linear Mechanics (ref. 5.5). 1 used this method on problems involving simple types of autopilots. With the eventual development of analog computers, such problems could be solved readily.
A final point that should be appreciated by the reader is the state of computational facilities at the time much of my work at Langley was conducted. Prior to about 1955, the only widely used computers were the slide rule and mechanical calculators such as the Marchant and Frieden that required several seconds to multiply or divide two numbers. To perform lengthy calculations, female employees called computers were employed to calculate results with these machines by following sheets that had the necessary steps listed in tabular form. Such calculation sheets are today called spreadsheets.
In making analytical studies of problems, a very desirable result was a closed-form solution, which means a formula in which numbers can be substituted for any particular case to determine the numerical value of the desired quantity. The determination of these closed-form solutions had been the objective of mathematicians and scientists for many years. In the aeronautical field, for example, Ludwig Prandtl determined the formulas for calculating induced drag of wings, Max Munk derived formulas for moments acting on ellipsoids in steady flow, and Theodore Theodorsen derived formulas for the lift and moments on oscillating wings. The formulas are considered to be in closed form if they give results in terms of known tabulated functions, such as trigonometric functions or Bessel functions. Most of these formulas and their originators became very famous because airplane designers could calculate results accurately for a configuration that at least approximated the one in which they were interested
With the advent of analog computers and later of high-speed digital computers, the need for closed-form solutions was reduced. Most calculations on the digital machines are made by numerical methods. One such method is the step-by-step method in which the response of a system to a disturbance or the trajectory of a vehicle is calculated a small increment at a time. Another numerical method is the Monte Carlo technique, in which many trial solutions are made to find the one with the best answer for the problem. Many problems of structures or fluid mechanics are solved by finite-element methods, in which the equations relating each small element of a large structure or flowfield are solved. These methods involve tremendous amounts of numerical computation, which modern electronic computers can handle in a very short time. Since I had none of these methods in my college education, my facility with such techniques is much less than that of students who grew up in the computer age.
With this background, I will describe some of the analytical studies that I conducted during my early employment at Langley.
Graphical Solution of the Quartic Equation
As pointed out in the introductory section, the solution for the motion of an airplane following a disturbance was very tedious because of the numerical calculations required. One of the problems encountered was solving for the roots of higher degree equations. A fourth-degree equation resulted in the solution for the longitudinal motion of a rigid airplane. For many problems, this equation could be reduced to a second-degree, or quadratic, equation by considering the airspeed constant, a valid assumption if short-period maneuvers were being considered. The lateral equations likewise result in a fourth-degree equation if certain simplifying assumptions are made. These fourth-degree equations, also known as quartics, are called the characteristic equations for the systems, and solving for the roots of these [41] equations is the first step in calculating the response of an airplane to controls.
Quadratic, cubic, and quartic equations may be written as follows:
High school students of mathematics are familiar with the formula for the solution of a quadratic equation. Formulas also exist for the solution of cubic and quartic equations, though they are considerably more complicated than those for a quadratic. For this reason, graphical methods or methods of successive approximations have been sought for the solution of these equations. In a paper on propeller governors, Herbert K. Weiss presented a set of graphs for the solution of cubic equations in terms of two parameters calculated from the coefficients of the equation (ref. 5.6). Shih-Nge Lin, in appendix I of his MIT thesis, described a method that he developed for a solution of quartic and other higher order equations (ref. 5.7). This procedure is a method of successive approximations that uses repeated long division. This method, though it was considerably quicker than the classical techniques, was still quite time-consuming without considerable practice. I therefore attempted to develop a set of charts for the quartic similar to those of Weiss for the cubic. A brief summary of this method and a sample copy of the charts are presented as appendix IV. These charts were never published. I used them to some extent in my work, but found that the time to solve for the parameters used in the charts and then to return the solution to an expression involving the original variables, required about as much time as Lin's method. This method is presented mainly because it was my only excursion into pure mathematics. With the development of high-speed computers, of course, the engineer no longer has to be concerned with these calculations because readily available computer programs can solve quartic or even much higher degree equations in a fraction of a second. As late as 1962, however, there was still interest in simplified methods for determining the roots of algebraic equations, as shown by a report by James W. Moore and Rufus Oldenburger of which I have an unpublished copy. This report presents a systematized procedure similar to Lin's method and analyses problems of convergence for cases, such as unstable roots, in which convergence of the method may be slow. Oldenburger was a well-known expert on servomechanisms before the subject of automatic control became a favorite subject for control theorists.
The airplane with fixed controls in the unstalled flight regime, fortunately, is beautifully linear. That is, all the aerodynamic forces and moments increase in proportion to the magnitude of the displacement or angular velocity, even down to very small magnitudes of motion. As a result, the measured motion of airplanes had been found to be closely predicted by the theory. Another favorable feature of linear systems is that a given solution is applicable to all magnitudes of motion. An increased initial disturbance simply increases all quantities involved in proportion to the disturbance without changing their time dependence. This condition no longer exists with nonlinear systems, and a separate solution has to be obtained for each magnitude of motion.
Solution of a Nonlinear Problem
The problem of snaking oscillations, discussed in the previous chapter, is an example of a nonlinear problem. The nonlinearity arises because of friction in the rudder control system. The rudder, instead of moving in proportion to the motion of the airplane, sticks until the aerodynamic hinge moments exceed the friction. Then the rudder starts to move under the influence of both the aerodynamic forces and the friction force.
[42] By the time the snaking motion had been explained, I was familiar with another method to analyze nonlinear systems. Robert T. Jones, while analyzing the stability of some of the first guided missiles developed by the Navy, had devised a technique assuming that the vehicle performs a sinusoidal oscillation. The actual control motion resulting from this sinusoidal motion was then calculated or measured experimentally. He then calculated the response of the system under the assumption that the work done on the vehicle per cycle must be the same as that done by the action of the control, and the angular impulse imparted to the vehicle over a half cycle must equal the change in angular momentum of the vehicle caused by the operation of the control. The work and momentum relations give two equations from which the frequency and amplitude of a constant-amplitude oscillation may be calculated. The expressions for the work and momentum imparted by the control may be shown to be related to the lowest order cosine and sine components of the motion of the control when expressed as a Fourier series. Jones' technique was therefore equivalent to the frequency-response method, which was then generally unknown to aeronautical engineers. It was later found that it had been developed to an extensive degree with different notation by electrical engineers for studying the stability of feedback amplifiers. I was intrigued by this method, because it appeared that the exact form of the control surface motion when subject to static friction would be important. As a result, I attempted to analyze the snaking problem by this method. I first tackled the problem of aileron snaking, a steady oscillation in roll caused by friction in the aileron control system. This problem had never been encountered in practice because of the much larger damping of the roll subsidence mode of an airplane when compared with the damping of the lateral oscillation. A much larger tendency for the controls to float against the relative wind would have been required to produce an aileron oscillation than to produce a rudder oscillation. The problem was easier to formulate, however, because of the simpler equations governing the rolling motion. Even in this case, calculation of the aileron motions subject to the control friction and the aerodynamic moments on the aileron proved to be quite difficult. After working on the problem as a part time activity for several months, I succeeded in obtaining a closed form solution. A typical example of possible steady oscillations for two cases is given in figure 5.1. In the first case, the friction is enough to cause the aileron to stick during part of the cycle; in the second case, the aileron motion is continuous.
After going through all this work, I concluded that it was hardly worthwhile to attempt a similar solution for the more practical case of the rudder snaking oscillation. The results on the aileron oscillations were never published. Many years later, however, a problem was encountered of control of a missile in which the ailerons were operated by a gyroscope sensing rolling velocity. This problem was exactly equivalent to the problem that I had solved. The only difference was that the hinge moments proportional to rolling velocity were applied to the ailerons by the gyroscope rather than by an aerodynamic floating tendency. It was possible to use my analysis to obtain some design data for the missile. As a result, my analysis was not completely wasted after all.
The frequency-response technique is based on the assumption that although the control motion is nonsinusoidal, the response of the controlled vehicle is very nearly sinusoidal because the relatively large inertia and slow response of the vehicle prevents it from responding appreciably to the higher frequencies in the irregular control motion. This assumption justifies the calculation of the control motion that is based on a sinusoidal vehicle motion. In most discussions of this method, however, the use of the lowest harmonic of the control motion in calculating the vehicle response is not justified on physical grounds. Jones' independent approach, which first based the response of the vehicle on the momentum and work relations over a....

[44] ....cycle and later showed that these relations gave expressions corresponding to the lowest order harmonic of the control motion, gives a clear physical interpretation of the application of the frequency-response method to nonlinear systems and provides insight into the accuracy of the method for this purpose. In the ensuing years, methods based on the frequency-response method were extensively used at the NACA and elsewhere for design of missile control systems, determination of stability derivatives from flight data (now called parameter identification), and calculation of response to turbulence. The many developments of the method, such as Bode plots, the Nyquist criterion, and root locus techniques gradually became part of the aeronautical engineer's mathematical equipment.
Computation of Lateral Oscillation Characteristics
As pointed out previously, the calculation of the airplane modes of motion prior to the introduction of high-speed computers was a tedious process. The solution for the lateral modes of motion requires the solution of a quartic equation. To simplify the process, approximate procedures were developed. Simple expressions could be found for the characteristics of the spiral mode and the roll subsidence, but accurate calculation of the Dutch roll usually required solution of the complete equations. A graphical method, first developed in 1937 by R. K. Mueller of MIT, proved to be applicable to the Dutch roll mode. Mueller first used the method in developing the world's first electric analog computer, which he used for the calculation of the longitudinal motion of airplanes (ref. 5.8). He discovered, however, that once he had perfected the graphical method, he no longer needed the analog computer. Later the method was discovered independently by K. H. Doetsch and W. J. G. Pinsker, two former German engineers working at the Royal Aircraft Establishment (RAE) after WW II. They called the technique the time vector method and applied it to calculation of the Dutch roll mode. Later, W. O. Breuhaus of the Cornell Aeronautical Laboratory, on a visit to England, learned the method from the RAE engineers and publicized it in the United States (ref. 5.9). Prof. E. E. Larrabee of MIT, on a sabbatical leave at the NACA Langley laboratory, showed how the method could be used to calculate stability derivatives from measured lateral oscillation data (ref. 5.1 0).
The time vector method is based on the fact that in a free oscillation of a linear system, the variables involved always maintain the same ratios of amplitude and the same relative phase angles. The variables can be depicted on a polar vector diagram in which the amplitude of each variable is shown by the length of an arrow and the relative phase angles by the directions of the arrows. In a damped oscillation, the diagram would rotate and shrink, but always maintain the same relative magnitude and angular separation of the vectors. All information about the motion could therefore be obtained by taking a snapshot of the diagram, which showed it at a given instant of time. Fixed relations can also be calculated between displacement, velocity, and acceleration of any quantity based on the frequency and damping of the motion.
The equations of motion involve a series of terms in each equation. In the lateral case, the equations involve roll, yaw, and sideslip. Each term in the equations consists of one of these variables multiplied by a stability derivative which gives the variation of the force or moment with the variable. For example, a typical term might be rolling velocity multiplied by the derivative, variation of rolling moment with rolling velocity. All the rolling moments contributed by the different variables are added up to give the total rolling moment, which must equal the inertia in roll multiplied by the rolling acceleration. This term may be placed on the opposite side of the equation, which results in a sum of terms equal to zero.
[45] Since the derivatives are constants, each term varies with time in exactly the same way as the variable in that term. A polar vector diagram may therefore be drawn that includes each term in the equation with the tail of each vector coinciding with the head of the previous one. Since the terms add up to zero, the diagram must close. A similar diagram may be drawn for each of the variables: roll, pitch, and yaw.
In using the method, the quantities to be determined are the frequency, the damping, and the ratios of the variables. Usually a value of unity is assigned to one of the variables, such as sideslip, and the ratios of roll to sideslip and yaw to sideslip are to be determined. All other quantities must be known beforehand or estimated. Usually the frequency is estimated first from the simple relation of the frequency of the airplane oscillating with a single degree of freedom in yaw. The damping is assumed to be zero. From these assumptions, initial values of the ratios of the variables may be determined by the closure of the vector diagrams. Then, with these values, more accurate values for the frequency and damping may be determined and these values used in a second attempt. The method is therefore an iterative procedure. Usually the convergence is very rapid and requires only two or three iterations to reach a solution.
The time vector method gained considerable recognition because of the rapidity of solution and because it gave a useful physical picture of the relations between the variables in an oscillation. It occurred to me, however, that the graphical work might be avoided if a similar convergent procedure could be performed analytically. After several trials of different methods, I discovered that the equations could be solved for the ratios of the variables. Therefore, a similar iterative procedure could be set up starting with an assumed value of frequency and a value of damping of zero, to calculate values of the ratios of the variables. I discussed this method with Bernard B. Klawans, an engineer in my branch. He worked out the details of the procedure. As it turned out, with the initial values of frequency and damping, the ratio of roll to yaw could be calculated. Then with this value and the values of frequency and damping, the ratio of sideslip to yaw could be calculated. Finally, with these two ratios, a quadratic equation for the Dutch roll root could be obtained that gave new values for the frequency and damping. Mr. Klawans checked the results for a wide range of variables. The calculations could be carried out readily on a slide rule or mechanical calculator, and the solutions converged very rapidly, usually within two or three iterations. Mr. Klawans published a Technical Note on the results (ref. 5. 11).
I thought this procedure was a worthwhile improvement over the graphical method, but evidently it arrived on the scene too late. In that period, about 1956, analog computers were already available, and some digital computers had been introduced. Soon programs were provided to solve the equations of motion in a matter of seconds. As a result, the approximate procedures that had occupied the efforts of stability and control engineers for many years fell into disuse.