VOL. 8, NO. 11, NOVEMBER 1971 J.

AIRCRAFT 885

An Approximate True Damping Solution of the Flutter

Equation by Determinant Iteration
HERMANN J. HASSIG*
Lockheed-California Company, Burbank, Calif.

The difference between a true damping, or rate-of-decay, solution of the flutter equation
and the structural-damping-type of solution is highlighted. True damping solutions are
possible if unsteady aerodynamics can be expressed in terms of the complex variable p. If
the aerodynamics are given at discrete values of the reduced frequency fc, an approximate de-
termination of the true damping is possible by assuming that the aerodynamic forces for har-
monic motion are a good approximation for the cases of slowly increasing or decreasing ampli-
tude. A determinant iteration method for obtaining the solution is presented. Results ob-
Downloaded by CARLETON UNIVERSITY LIBRARY on June 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/3.44311

tained by different methods of solving the flutter equation are compared.

Nomenclature p, pQ = air density, air density at sea level

a- = air density ratio p/po
quasi-steady aerodynamic matrices w = circular frequency

[#2] , [.Bi] , [.Bo] = aerodynamic matrices denned by Eq. (8) The p and k Methods of Solution
(A(ik)],[
[Aa(k)] complex aerodynamic matrix numerically of the Flutter Equation
available at any discrete value of k
[A(p)] aerodynamic matrix that is an explicit func- HE flutter equation can be written in the general
tion of p form
[D] viscous damping matrix
matrices relating the control system to the = 0 (1)
[DA] , [Ds] ,
degrees of freedom [K] defines the elastic characteristics by relating generalized
[De],[Dz] matrices expressing local slope and local ver-
forces {Q} to generalized displacements {q}
tical displacement in terms of {q}
\[F(p,k)]\ flutter determinant = [Q\ (2)
[K] stiffness matrix
[M] inertia matrix [If ] defines the inertia characteristics by relating the inertia
[M] inertia matrix including effective apparent forces {Qin} to the generalized accelerations
inertia of air
{Q} column of forces in} = -lM}{q} (3)
{ q} column of degrees of freedom
column of angles of attack [A (p) ] defines the unsteady aerodynamic forces through the
constants defining the lag function equation
an,an+i successive amplitudes of oscillatory motion {Qaero} = ipV*[A(p)]{q] (4)
c reference chord
F,G scalars In Eqs. (1) and (4) p is the nondimensional differential opera-
/ frequency in cps tor (c/V)(d/dt). If the aerodynamic forces can be expressed
scalar as a sufficiently simple function of p, Eq. (1) defines a poly-
structural damping nomial in p with real coefficients. For nonzero solutions for
gravitational acceleration q the determinant formed by the matrix coefficients in Eq. (1)
transfer functions for hydraulic controls and must be equal to zero. For a given value of the speed V, the
automatic control systems; explicit func- determinant can then be solved directly for p. This leads to
tions of s
k reduced frequency coc/F conjugate complex roots
p differential operator (c/V)(d/dt) ik (5)
p = 5 + ik iterated values for the root of the flutter equa-
tion k defining a nondimensional reduced frequency coc/7 and 7
nondimensional time Vt/c defining a rate of decay
V true airspeed
7 true damping coefficient 7 = (l/2r) ln(an+1/«») (6)
X = \r i\i = complex root of flutter equation
where an and an+\ are the amplitudes of successive cycles.
This method of solving the flutter equation will be called
Received February 8, 1971; revision received July 6, 1971. the p method. It can be used in the case of quasisteady aero-
The author wishes to acknowledge Robert B. Neveceral of the dynamics when [A (p) ] assumes the simple form
Scientific Computer Services Division of the Lockheed-California
Company who did all the computer programing for this paper. [Ai]p+ [A0] (7)
The work presented herein was conducted as part of Lockheed-
California Company independent research and development. but also in the case of simple forms of unsteady aerody-
Index category: Aeroelasticity and Hydroelasticity. namics.
* Research and Development Engineer, Commercial Engineer- Mazelsky and O'Connell1 formulated approximate un-
ing. Member AIAA. steady aerodynamics that accounted for apparent mass, aero-
 H. J. HASSIG J. AIRCRAFT
COMPARISON BETWEEN p-METHOD & k-METHOD
FOR A TWIN-JET TRANSPORT AIRPLANE
by formulating the flutter equation for a twin-jet transport
k = 2.0 . 1.0
according to Eq. (8) and solving it according to the p method
and the k method.
The upper part of Fig. 1 shows the frequency-vs-speed rela-
tion obtained with the two methods. Each of the curves
numbered 1-4 corresponds to a root of the flutter equation
and an associated flutter mode. The lower part of Fig. 1
shows how the damping obtained with the two methods varies
with speed.
For the k method g is plotted; for the p method the com-
parable quantity
o 2y = (1/vr) ln(an+1/an) (Ref. 6)
-0.1 The numbering of the curves corresponds to the numbering in
the upper part of the figure. One can see that the p method
-0.2 solution shows flutter in the first mode; the k method solution
shows flutter in the second mode. Even the modal coupling
-0.3
300 400
seems different. Yet at the flutter speed, both methods of
Downloaded by CARLETON UNIVERSITY LIBRARY on June 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/3.44311

V -KEAS solution give identical results.

Fig. 1 Comparison between p method and k method for It is generally conceded that it is desirable to formulate and
a twin-jet transport airplane. solve the flutter equation such that the solution leads to a
value for the rate of decay. Ideally, this requires the formu-
lation of the unsteady aerodynamics matrix as a function of
dynamic lag, and aerodynamic coupling between different the complex variable p = yk + ik. Some of the approximate
strips on the wing and which were a relatively simple function formulations of the aerodynamic matrices have made this
of p. With their formulation, the flutter equation takes the possible.1'4-5 An attempt to develop Theodorsen-type aero-
form dynamics as a function of p = yk -f ik led to a dispute as to
+ p[Al] + [A,]] + whether the formulation was valid for a motion with decaying
2]
amplitude.7"9 The only published application of generalized
+ [B^]{q] = 0 (8) Theodorsen aerodynamics is by Stiimke.10 No such attempt
The quantities bi and &2 are from the two-exponential lag for kernel function aerodynamics or supersonic Mach box
function aerodynamics has come to the attention of the author.
In general, then, when one wants to work with exact
1 - aie~b^s + atf-b*s (9) theoretical aerodynamics one must work with a formulation
for harmonic motion and devise approximate methods to de-
The matrices A 0 ,Ai,A 2 and B0)Bi}Bz contain the lift curve
termine the rate of decay.
slope, aerodynamic influence coefficients, the aerodynamic Zisfein, Frueh, and Miller11'12 have shown that under sim-
center and rotation center (corresponding to the f chord plifying assumptions, the rate of decay as a function of speed,
point in two-dimensional theory) for each strip, and the damp- assuming zero structural damping, can be obtained from the
ing in pitch. Equation (8) has been used at the Lockheed- traditional /c-method solution. If g, co, and V are found with
California Company since the middle fifties and has been the k method, then, according to Ref. 12
routinely solved by determinant iteration since 1960 for roots
p = jk + ik of interest.!
More sophisticated formulations of the aerodynamics, such
as follow from the kernel function or doublet lattice approach
or the supersonic Mach box, lead to aerodynamic matrices Landahl13 plots the real and imaginary part of the flutter de-
valid only for harmonic motion, p = ik. In that case, Eq. terminant as a function of co and obtains approximations for
(1) takes the form 7 and co. Such a plot must be made for each speed to obtain
complete frequency damping-velocity diagrams. Natke 2 also
- 4 [M]k* [K] - = 0 (10) describes a method of obtaining rate-of-decay from what is in
principle a /c-method type solution. These methods are ap-
This is the traditional American form of the flutter equation. proximate and indirect. Indirect in the sense that they re-
At chosen values of /c, complex roots for 1/T2, X r + i\i, are quire the availability of the /c-method solution.
found and interpreted as
The p-k Method

where g is the structural damping required for harmonic mo- Equation (1) suggests an approximate method of finding a
tion. This method of solving the flutter equation is here rate-of-decay type solution directly. Writing Eq. (1)_ in a
called the k method. form indicating:that the aerodynamic matrix is only available
Several authors have discussed the differences between the for harmonic motion, and computing [A(ik)]ior an estimated
p method and the k method and variations thereof.2"5 The value of k :
significance of the difference between the rate-of-decay [K] - (13)
damping, 7, found with the p method and the structural
damping, g} found with the k method is, in general, well under- one can solve for p = yki + iki, compute [A(iki)], solve Eq.
stood. The numerical differences have been demonstrated (13) again which leads to p = ykz + ikz, etc., until the imagi-
by direct calculation by Richardson, 4 Rodden and Stahl.5 A nary part of the solution equals the k value of the aerody-
result obtained by this author is shown in Fig. 1. It is found namics. This method of solution is here called the p-k
method. The rationale for this approach is that f or ^ sinu-
soidal motions with slowly increasing or decreasing amplitude,
t Credit for the original determinant iteration program goes
to R. F. O'Connell and G. E. Smith of the Lockheed-California aerodynamics based on constant amplitude are a good ap-
Company. proximation.
 NOVEMBER 1971 DAMPING SOLUTION OF FLUTTER 887
COMPARISON BETWEEN p-METHOD & p-k-METHOD
FOR A TWIN-JET TRANSPORT AIRPLANE
The method is shown in an actual flutter analysis by Irwin SOLID LINES - p-METHOD
and Guyett14! and is mentioned by Natke2 -and Dat and °g]- p-k-METHOD
Muerzec.3 Recently, Jocelyn Lawrence and Jackson15 de-
voted an ARC report to a comparison of methods of solving
the flutter equation and they presented the method of Ref.
14 in detail, calling it the "British Method" in contrast to the
"American Method" (k method). §
In Refs. 14 and 15, a graphical method is presented to
match the imaginary part of p with the k value of the aero-
dynamics. This author has shown the feasibility of finding
roots directly. Figure 2 compares p-k method results and p
method results obtained by applying Eq. (8) to the same case
of a twin-jet transport as in Fig. 1. Corresponding curves in
the upper and lower halves of the figure have like symbols.
At all rates-of-decay of significance the p-k method gives a
very good approximation of the true rate-of-decay as found
with the p method.
In judging the significance of Fig. 2 it must be considered 100 200 300 400 500 600
that the aerodynamics implied by Eq. (8), although contain-
Downloaded by CARLETON UNIVERSITY LIBRARY on June 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/3.44311

V - KEAS
ing most important aspects of unsteady aerodynamics, is in a Fig. 2 Comparison between p method and p-k method
rather simple form. It does include aerodynamic lag and for a twin-jet transport airplane.
steady-state aerodynamic coupling between strips, but not
the signal delay between strips. However, Fig. 2 represents
the most general comparison known to this author at the Eq. (14) the values
moment.
In view of the good agreement shown in Fig. 2, a computer Fi=\[F(pM]\ Fi = \[F(PM (17)
program was developed that makes it possible to find roots
according to the p-k method by determinant iteration. are determined. The Reglua Falsi method gives a first
iterated value for p :
Matrix Iteration Applied to the p-k Method p3 = (pzFi - piFz)/(Fi - F2) (18)
When the p-k method symbolized by Eq. (13) is pro-
gramed for the digital computer, significant advantages are The process is repeated according to the recurrence formula
gained by generalizing the equation. The equation used in
the Lockheed-California Company program is F,+1) (19)
I F 2 until a specified degree of convergence is attained. From the
- [M] -2 p2 + [D] - p + (1 + ig)[K] -
go c c converged root pc — 5C -+- ikc, the frequency and damping
can be computed
V \
- HA [DA] - HB [DB] -
~cP) Vkc
(20)
{q} = 0 (14)
To complete one frequency-damping-velocity diagram for
or several modes, one initial trial for each mode, at the first
speed only, must be input. This trial, pz, may be given as
[F(p,K)]{q] =0 (15) P2 = ^2 + ik<> if available from an earlier computation, or as
[D] represents viscous damping; e.g., due to flutter dampers.
HA, HB, and Hc are transfer functions lor hydraulic controls p2 = 0 + (2irfc/V)i (21)
and automatic control systems; they are related to the de-
grees of freedom by [DA], [DB], and [Dc]^ I/go and/L are where/is the natural vibration frequency of the mode sought.
scalars, introduced for convenience. When using Eq. (21), the first speed, 7, should be sufficiently
All matrices in Eq. (14) are real and uniquely defined, ex- small such that Eq. (21) is a reasonably good trial for p2.
cept [A(A;)], which is complex and must be given for a suffi- The program computes pi as
cient number of k values. Equation (14) is solved at several
values of V and a, or combinations thereof, for complex roots pi - -Fk* + (22)
p associated with modes of interest. Modes of interest are where usually F = 0.01 and G = 1.00.
determined from vibration analysis or from previous flutter When all the desired roots at the first speed are found, they
analyses. are used in determining the initial trials for the next speed, V
The process of determinant iteration is completed mode by
mode for one speed and then at successive preselected speeds. Pi = (V/V)p* P2 = (V/V)Pc (23)
For one mode at one particular speed, the process is started by
initial trials for p: and repeated application of Eq. (19) leads to converged roots
at the second speed
Pi = <5i + iki p2 = 52 + ikz (16)
] and [A(fe)] are computed by interpolation. Using PC = dc + ike (24)

t This author's attention was directed to this method by The trials for the third speed, V, are
P. R. Guyett when he was visiting the Lockheed-California Co.
§ This author favors identifying methods by names that say Pi = (V/V)pe p2 = (V/V)pc (25)
something about the method, if it can be done without using un-
duly long names. and similarly for the following speeds.
 H. J. HASSIG J. AIRCRAFT
COMPARISON BETWEEN p-k-METHOD & k-METHOD
FOR A HORIZONTAL STABILIZER WITH ELEVATOR
D / vs V and 2y vs V plots can be obtained from the CALCOMP
digital incremental plotter.

Discussion
A good deal of experience with the Determinant Iteration
p-k Method Program has been gathered at the Lockheed-
California Company. Figure 3 shows one of the early check
cases for which an earlier /c-method solution was available.
The /c-method solution is represented by the solid lines with
computed points marked. Corresponding curves in upper
and lower half of the figure are numbered 1-4. The p-k
method solution is represented by the symbols and where it is
significantly different from the /c-method solution, by dashed
lines. A case of a horizontal stabilizer with elevator is chosen
for which the original k solution can be easily interpreted in-
correctly even though the k values at which the flutter equa-
tion is solved are reasonably close together. A sufficient
Downloaded by CARLETON UNIVERSITY LIBRARY on June 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/3.44311

number of check cases has been done to show the agreement

at the flutter speed between k method and p-k method.
It is thought that the p-k method is economically competi-
Fig. 3 Comparison between p-k method and k method
for a horizontal stabilizer with elevator.
tive with the k method, although no precise comparisons have
been made. Much depends on details of the computer system
available. The p-k method has two distinct advantages,
Features of the Lockheed-California however, the frequency-damping-speed plots are more easily
Company Program interpreted and frequency-damping-speed plots can be ob-
tained in which the solution at each speed is for the density in
The program, Determinant Iteration p-k Method, was pro- the atmosphere corresponding to that speed and the Mach
gramed by R. B. Neveceral for the IBM 360 Model 91 as number. Thus, computer runs for several different constant
part of the FAMAS system.16 The existing version is de- densities are avoided.
signed such that each iteration can be done largely within core. In developing the program, elements of the [A(/c)] matrix
Therefore, the following restrictions have been imposed: a) were plotted vs /c, and surprisingly irregular curves were found
maximum order: 50 X 50; b) use four [A(/c)] matrices for for the case in which a 64-degree-of-freedom system was re-
the interpolation of the aerodynamic matrices; c) the duced to a modalized 15-degree-of-freedom system, using
matrices [DA], [Ds], and [Dc] have a maximum of 20 nonzero natural vibration modes. This has given rise to the hunting
elements; and d) HA, HB, and He each are the quotient of 8th mentioned earlier. When eliminating hunting, by retaining
degree polynomials in (V/c)p. the same four aerodynamic matrices as long as the imaginary
The four aerodynamic matrices used for each interpolation part of p does not wander away too far, the problem of non-
are determined by the value of the imaginary part of p used uniqueness of the aerodynamics is introduced. A compro-
to evaluate | [F(p,k) ] \ . During the first three iteration steps, mise can be found in which hunting is virtually avoided, and
the aerodynamic matrices are chosen such that two of their k nonuniqueness of the aerodynamics is restricted to small
values are above and two are below the imaginary part of p. ranges of k values (say, 5% of each interval at each end of the
From then on, the same four aerodynamic matrices are used; interval).
even if during the iteration process the imaginary part of p Possibly better behaved aerodynamic matrices are found if,
wanders outside the range between the inner two k values, as instead of writing Eq. (4), the aerodynamic forces are ex-
long as it does not wander too far into the adjacent range.
pressed as
This feature was included after it was found that, when the
final k value of the solution is close to one of the input values, {Qaero} = ±PV*[Aa(k)]{ a} (26)
"hunting" may occur. Namely, going from one iterated
value to the next, the imaginary part may demand a different from which
set of aerodynamic matrices and the next iterated value may (27)
call the earlier four aerodynamic matrices back.
For the determinant iteration to work satisfactorily, the However, this increases the number of matrices that must be
initial trials that start the iteration must be reasonably good. kept in store for the iteration; as a result, the maximum order
Thus, when going from one preselected speed to the next, a that can be handled will be less than 50.
root may vary so much that the trials provided by the pro- At low speed and high frequency, the k values are high and
gram do not lead to convergence. In that case, for the mode often lead to extrapolation for [A (k) ], which in turn may lead
concerned, the program cuts the interval to the previous to nonconvergence or convergence on the wrong roots. To
speed in smaller subintervals. avoid problems, an aerodynamic matrix for an arbitrarily high
As a result of the previous feature, speed intervals may be k value is included (say, k = 50 or 100). The [A(k) ] for that
chosen rather large, such that when flutter occurs there are value is derived from piston theory or made equal to the
not enough points to determine the flutter speed accurately regularly computed [A (k) ] matrix at the highest k value.
by fairing a curve. To get a better definition of the flutter The present program makes it possible to find control sur-
speed, the program will solve Eq. (14) at additional speeds in face rotation modes that in the k method often escape atten-
the flutter region any time a mode becomes unstable. tion because, in the frequency-vs-speed diagram, they fall
There is an option to run the program at constant a (a between two constant k lines. Also, the stability modes,
maximum of five a values per case) or at cr varying with except the phugoid, can be found at the same time that the
speed. If in the latter case a varies as defined by a constant flutter modes are found.
Mach number in the atmosphere, the program computes the The program can be expanded to include computing the
a values corresponding to the velocities chosen. The a values effect of parameter variations on the frequency and the damp-
for any other speed-vs-0- relationship must be listed on a con- ing at one speed. Once the frequency and damping at that
trol card. The output of the program is arranged such that speed have been determined, say for [Mi], [Di] or [Ki] the
 NOVEMBER 1971 DAMPING SOLUTION OF FLUTTER EQUATION

roots can be used as first trials to start the iteration at the the So-Called 'Reduced Frequency Scanning' Method," La
same speed for a slightly different matrix [M2], [D%] or [jK 2 ]« Recherche Aerospatiale, No. 133, Nov.-Dec. 1969, pp. 41-43.
4
This may have advantages in the case of the flutter analysis Richardson, J. R., "A More Realistic Method for Routine
Flutter Calculations," AIAA Symposium on Structural Dynamics
of airplanes with a large number of possible external wing store
Downloaded by CARLETON UNIVERSITY LIBRARY on June 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/3.44311

and Aeroelasticity, Aug. 30-Sept. 1, 1965, Boston, Mass.

combinations, or for other parameter studies. The feasibility 5
Rodden, W. P. and Stahl, B., "A Strip Method for Prediction
of this process was recently demonstrated for the variation of of Damping in Subsonic Wind Tunnel and Flight Flutter Tests,"
the viscous damping [D]. Journal of Aircraft, Vol. 6, No. 1, Jan.-Feb. 1969, pp. 9-17.
6
Determinant iteration was used for finding the flutter roots Scanlan, R. H. and Rosenbaum, R., "Further Consideration
because it was readily available. It can be used without first of Linear Damping," Introduction to the Study of Aircraft Vibra-
putting the flutter equation in the canonical form and the tion and Flutter, 1st ed., Chap. Ill, Macmillan, New York, 1951,
flutter determinant can be an intricate function of p. (Thus, pp. 85-88.
7
determinant iteration can be used if aerodynamic approxima- Luke, Y. L. and Dengler, M. A., "Tables of the Theodorsen
Circulation Function for Generalized Motion," Journal of the
tions are complicated formulas of p, such as those proposed Aeronautical Sciences, Vol. 18, No. 7, July 1951, pp. 478-483.
by Richardson.4) However, other methods of finding the 8
Van de Vooren, A. L, "Generalization of the Theodorsen
roots may be modified to suit the needs of the p-k method, and Function to Stable Oscillations," Journal of the Aeronautical Sci-
it has already been suggested17 that a modified power method ences, Vol. 19, No. 3, March 1952, pp. 209-211.
9
(matrix iteration) may be used successfully. Laitone, E. V., "Theodorsen's Circulation Function for
Generalized Motion," Journal of the Aeronautical Sciences, Vol.
19, No. 3, March 1952, pp. 211-213.
Conclusions 10
Stiimke, H., "Frequenzgang des ebenenen starren Trag-
fliigels in inkompressibler Potentialstromung und eines einf achen
1) An approximate but direct method of finding rate-of-
Ersatzschwingers," Zeitschrift fur Flugwissenschaften, Vol. 17,
decay type solutions of the flutter equation, even when the No. 12, Dec. 1969, pp. 425-433.
aerodynamic forces are available for harmonic motion only, 11
Zisfein, M. and Frueh, F. J., "A Study of Velocity-Fre-
is presented. 2) A determinant iteration method is used for quency-Damping Relationships for Wing and Panel Binary
solving the flutter equation which is generalized to include Systems in High Supersonic Flow," Rept. 59-969, Oct. 1959, Air
viscous damping, structural damping and several transfer Force Office of Scientific Research, Arlington, Va.
12
functions representing automatic control systems. 3) The Frueh, F. J. and Miller, J. M., "Prediction of Dynamic Re-
validity of the method is demonstrated by comparing solu- sponse from Flutter Analysis Solutions," Rept. 65-0952, June
tions obtained by different methods. 4) Matrix iteration 1965, Air Force Office of Scientific Research, Arlington, Va.
13
Landahl, M. T., "Graphical Technique for Analyzing Mar-
(power method) and other methods, may, after modification,
ginally Stable Dynamic Systems," Journal of Aircraft, Vol. 1, No.
provide alternate means of finding the flutter roots according 5, Sept.-Oct. 1964, pp. 293-299.
to the p-k method. 14
Irwin, C. A. K. and Guyett, P. R., "The Subcritical Re-
sponse and Flutter of a Swept Wing Model," Tech. Rept. 65186,
References Aug. 1965, Royal Aircraft Establishment, Farnborough, U.K.
15
Jocelyn Lawrence, A. and Jackson, P., "Comparison of Dif-
1
Mazelsky, B. and O'Connell, R. F., "Transient Aerodynamic ferent Methods of Assessing the Free Oscillatory Characteristics
Properties of Wings: Review and Suggested Electrical Repre- of Aeroelastic Systems," Current Papers, C.P.N. 1084, Dec. 1968,
sentation for Analog Computers," LR 11577, July 1956, Lockheed Aeronautical Research Council, London, England.
16
Aircraft Corp., California Div., Burbank, Calif. "Flutter and Matrix Algebra System—F AM AS—Manual,"
2
Natke, H. G., "Bemerkungen zu der angenaherten Losung LR 23657, April 20, 1970, Lockheed-California Co., Burbank,
des klassischen Flatterproblems," Zeitschriftfiir Flugwissenschaf- Calif.
17
ten, Vol. 15, No. 11, Nov. 1967, pp. 425-431. MacNeal, R. H., private communication, Dec. 17, 1970,
3
Dat, R. and Meurzec, J. L., "On the Flutter Calculations by The MacNeal-Schwendler Corp., Los Angeles, Calif.