Integration Methods and Closed Forms

Digital computer programs for airplane stability and control time history calcula­tions perform step-by-step integration of the equations of motion. The usual form of the complete equations of motion for numerical integration on a digital computer is 12 simul­taneous nonlinear first-order differential equations. Three of the equations produce linear position coordinates, or state components, three produce attitude angles (if Euler angles are used), three produce linear velocity components, and three produce angular velocity components. The 12 airplane coordinates of motion are referred to as the airplane’s state vector.

Accurate, efficient integrating algorithms were a subject of interest among applied math­ematicians for many years before stability and control engineers needed them for computer programming. A well-known text that compares the properties of many integrating algo­rithms is Introduction to Numerical Analysis by F. B. Hildebrand, published by McGraw-Hill in 1956.

A fair generalization is that choice of an integrating algorithm is a trade-off between simplicity, which affects calculation speed, and accuracy. The simplest algorithms, such as Eulerian or “boxcar” integration, require just one calculation pass per computing interval, but they accumulate systematic errors as the time history calculation goes forward. In Eulerian integration, a coordinate such as pitching velocity is projected forward to the next time interval simply by adding to the present value the product of the present value of pitching acceleration and the time interval length, which is usually of the order of 0.05 second. In general terms, the state vector at the next time interval is the current state vector plus the product of the state derivative vector and the time interval.

More accurate integration requires the calculation of intermediate values in order to take the same time step, adding to the computing time but improving accuracy. The best known accurate integration method, in effect a standard for stability and control time his­tory calculations, is the fourth-order Runge-Kutta method. This method can be adapted in FORTRAN to the integration of multiple states, such as the 12 airplane coordinates or state vectors (Melsa and Jones, 1973, and Figure 18.12).

While digital computers became available in engineering offices for stability and control time history calculations around the time of the inertial coupling crisis, or 1950, it was not until many years later that computing speed had increased to the point that the calculations could be made in real time and could thus support flight simulation. One of the earliest such applications was at the Ling-Temco-Vought plant in Arlington, Texas, in the late 1960s. An all-digital flight simulator that went on-line a little later was Northrop Aircraft’s Large Am­plitude Simulator, which progressed from analog to hybrid to all-digital in late 1975. With the introduction of real-time digital flight simulation accurate, but slow, integration methods such as the fourth-order Runge-Kutta routine have become something of an obstacle. There is a premium on the development of fast integration methods that still have a fair degree of accuracy Fast but accurate integration methods have been developed all over the United States to meet that need: methods generally starting with a classical scheme and modified by mathematical tinkerers.

The Adam-Bashford method was the starting point for the algorithm used for the projector gimbals in the Northrop Large Amplitude Simulator. A different set integrates the airplane equations of motion. Another integration method developed specifically for flight simulation modifies the second-order Runge-Kutta method, replacing the second state derivative vector calculation with a prediction based on a weighted average of previous mid – and endframe values. The modified second-order Runge-Kutta method seems to be almost as accurate as

SUBROUTINE INTG1

C FOURTH-ORDER RUNGE-KUTTA INTEGRATION DIMENSION Y1 (40),E1 (40),E2(40),E3(40) TSTEP = DT/2.

DO 2 l = 1,N 2 Y1 (I) =X(I)

CALL DERV1 DO 4 l= 1,N E1(I)=TSTEP*F(I)

4 X(I) = Y1(I) + E1(I)

T = T + TSTEP CALL DERV1

DO 5 1 = 1,N E2(l) = DT*F(I)

5 X(I) = Y1(I) + .5*E2(I)

CALL DERV1

DO 7 l = 1,N E3(I) = DT*F(I)

7 X(l) = Y1 (I) + E3(l)

T = T + TSTEP CALL DERV1

DO 8 1=1,N

8 X(l)=Y1(l) + (E1(l) + E2(l) + E3(l)+TSTEP*F(l))/3 RETURN

END

Figure 18.12 FORTRAN digital computer subroutine for the integration of a state derivative vector x. This is the widely used fourth-order Runge-Kutta method. COMMON input-output statements have been removed for generality. (From ACA Systems, Inc. FLIGHT program)

the fourth-order Runge Kutta, while requiring only one calculation of the state derivative vector per interval.

Aerodynamic forces and moments are involved in the calculation of the state derivative vector. This requires huge amounts of table lookup on modern flight simulations that cover large Mach number, altitude, and control surface position ranges and uses computer time more than any other part of the computation. Thus, a single calculation of the state derivative vector, as in the modified second-order Runge-Kutta method, is very efficient for real-time digital flight simulation. A modified second-order Runge-Kutta method was developed in 1972 by Albert F. Myers of NASA; it was then improved by him in 1978 for the HIMAT vehicle flight simulation (Figure 18.13).

Another important advance in digital flight simulation is the use of closed-form solutions for the first – and second-order linear differential equations that typically represent analog flight control elements, such as control surface valves and actuators. Closed-form solutions for these elements remove them from the state vector that has to be integrated, reducing the order of that vector to perhaps no more than is required by the airplane equations of motion themselves, or 12. The nonlinearities of control position and velocity limiting are easily represented. This technique is attributed to Juri Kalviste, although there may be other claimants to priority.