June 2015 

NASA/TP–2015-218765 

Re-Computation of Numerical Results 

Contained in NACA Report No. 496 

Boyd Perry, III

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NTRS Registered and its public interface, the 

NASA Technical Reports Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in both 

non-NASA channels and by NASA in the NASA STI 

Report Series, which includes the following report 

types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase of 

research that present the results of NASA 

Programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations. 

• TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 

preliminary or of specialized interest, 

e.g., quick release reports, working 

papers, and bibliographies that contain minimal 

annotation. Does not contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or 

co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI program, 

see the following: 

• Access the NASA STI program home page at 

http://www.sti.nasa.gov 

• E-mail your question to help@sti.nasa.gov 

• Phone the NASA STI Information Desk at 

757-864-9658 

• Write to: 

NASA STI Information Desk 

Mail Stop 148 

NASA Langley Research Center 

Hampton, VA 23681-2199 
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

 June 2015 

NASA/TP–2015-218765 

Re-Computation of Numerical Results 

Contained in NACA Report No. 496 

Boyd Perry, III

Langley Research Center, Hampton, Virginia 
Available from: 

NASA STI Program / Mail Stop 148 

NASA Langley Research Center 

Hampton, VA 23681-2199 

Fax: 757-864-6500 

Acknowledgments 

The author wishes to acknowledge the valuable contributions to this paper made by Mr. 

Thomas G. Ivanco of the NASA Langley, Aeroelasticity Branch and Mr. G. Lee Pollard 

of NCI Information Systems, Inc. 

Mr. Ivanco is a very talented engineer. While on a one-year active-duty assignment with 

the U.S. Army Reserve, during his off-duty time, he very generously performed all of 

the independent flutter calculations used as check cases in this paper. 

Mr. Pollard is a very talented graphic designer. He assembled the re-computed results 

generated by the present author, electronically scanned all the original figures from 

NACA 496, merged the two in figures 5 through 17, and produced all the final figures in 

this paper. 

Mr. Ivanco’s and Mr. Pollard’s skills and enthusiasm made my interactions with them not 

only very productive, but also very enjoyable. Many thanks to you both! 

The use of trademarks or names of manufacturers in this report is for accu rate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics and 

Space Administration. 
1

ABSTRACT

An extensive examination of NACA Report No. 496 (NACA 496), “General Theory of Aerodynamic

Instability and the Mechanism of Flutter,” by Theodore Theodorsen, is described.  The examination

included checking equations and solution methods and re‐computing interim quantities and all

numerical examples in NACA 496.  The checks revealed that NACA 496 contains computational shortcuts

(time‐ and effort‐saving devices for engineers of the time) and clever artifices (employed in its solution

methods), but, unfortunately, also contains numerous tripping points (aspects of NACA 496 that have

the potential to cause confusion) and some errors.  The re‐computations were performed employing the

methods and procedures described in NACA 496, but using modern computational tools.  With some

exceptions, the magnitudes and trends of the original results were in fair‐to‐very‐good agreement with

the re‐computed results.  The exceptions included what are speculated to be computational errors in

the original in some instances and transcription errors in the original in others.  Independent flutter

calculations were performed and, in all cases, including those where the original and re‐computed

results differed significantly, were in excellent agreement with the re‐computed results.  Appendix A

contains NACA 496; Appendix B contains a Matlab ® program that performs the re‐computation of

results; Appendix C presents three alternate solution methods, with examples, for the two‐degree‐of‐

freedom solution method of NACA 496; Appendix D contains the three‐degree‐of‐freedom solution

method (outlined in NACA 496 but never implemented), with examples.

I.  INTRODUCTION

In a year 2000 engineering note, Zeiler [1] pointed out that several of the foundational papers and texts

in aeroelastic flutter [2 through 6] contain numerical errors in some of their example problems.  It is not

surprising that such errors exist because – especially in the cases of references 2, 3, and 4, written in the

1930’s and 40’s – all calculations were performed “by hand” with pencil and paper, slide rules, and, after

they were invented, large table‐top four‐function electro‐mechanical calculators.

Because these foundational papers and texts are often used in graduate courses on aeroelasticity, Zeiler

recommended that an effort be undertaken to employ the computing resources available today to re‐

calculate the example problems and to publish the results so as to provide a complete and error‐free set

of example problems.  Following Zeiler’s recommendation, the purpose of this paper is to provide re‐

computed results for the example problems presented in reference 2:  NACA Report No. 496, “General

Theory of Aerodynamic Instability and the Mechanism of Flutter,” by Theodore Theodorsen (referred to

hereinafter simply as “NACA 496”).  Future papers will provide re‐computed results for references 3 and

4.

In order to perform the re‐computation of its numerical examples, an extensive study of NACA 496 was

undertaken by the present author.  The study entailed many parts, including reading and re‐reading the

paper, checking and re‐checking its equations, and checking and re‐checking its solution method.  In the

process, the present author discovered in NACA 496 computational shortcuts (time‐ and effort‐saving

devices for engineers of the time) and clever artifices employed in its solution methods.  Unfortunately,
2

the present author also discovered, via tripping, several tripping points (aspects of NACA 496 that have

the potential to cause confusion) and errors.  These items are discussed in the present paper.

In reference 1, Zeiler stated “One does not set about lightly to correct the masters ...”  Embracing this

notion, the present results have been carefully checked and re‐checked, and, as a further check, for

some of the numerical examples, independent calculations were performed by a colleague.  In all cases,

including those where the original and present results differed significantly, the independent

calculations were in excellent agreement with the present results, providing additional confidence in the

present results.

The arrangement of the remainder of the present paper is as follows:

Section II contains preliminaries, which are intended to aid in the understanding of the

remainder of this paper.   

Section III contains an overview of NACA 496, including its theoretical development and solution

method.   

Section IV describes the checks performed by the present author on the equations and solution

method of NACA 496.

Section V presents the re‐computation of quantities contained in Tables I‐IV of NACA 496.

Section VI presents the re‐computation of the example problems contained in Appendix II of

NACA 496.  

Section VII contains concluding remarks.

Appendix A contains a reproduction of NACA 496.  (Throughout the present paper the reader

will be referred to Appendix A frequently.)

Appendix B contains a Matlab® program that implements the two‐degree‐of‐freedom flutter

equations and the two‐degree‐of‐freedom solution method found in Appendix I of NACA 496.

Appendix C contains descriptions of, and example problems for, three alternate solution

methods for the two‐degree‐of‐freedom flutter equations presented in NACA 496.  

Appendix D contains a description of, and example problems for, the three‐degree‐of‐freedom

flutter equations presented in NACA 496.  
3

II.  PRELIMINARIES

As stated in the Introduction of the present paper, the information contained in this section is intended

to aid in the understanding of the remainder of this paper.   

Terminology

The present paper refers to sections and appendices within itself and sections and appendices within

NACA 496.  The present paper also refers to the authors of both papers.  To avoid confusion in regard to

which paper or author is being referred to, the following (unfortunately but necessarily) awkward‐

sounding convention is adopted:  for the present paper, the terms “present paper,” “this paper,” “main

body,” “Appendix A, B, C, or D of the present paper,” and “present author” are used; for NACA 496, the

terms “NACA 496,” “main text,” “Appendix I, II, or III of NACA 496,” and “author of NACA 496,” are used.

Symbology

The present paper retains the symbology of NACA 496.  Appendix A of the present paper (the

reproduction of NACA 496) contains a symbol list (pp. 9 and 10 of NACA 496).  Symbols not contained in

the list are defined as they are used in NACA 496 and as they are used in the present paper.

Tripping Points and Errors

As used in the present paper, the term “tripping point” is defined to be an aspect of NACA 496 that has

the potential to cause confusion for a reader of NACA 496.  Tripping points that the present author

encountered while reading NACA 496 are identified throughout the present paper and are collected in

Table I of the present paper.

Errors in some of the equations, tables, and figures of NACA 496 were found by the present author.  

These errors are identified throughout the present paper and are collected in Table II of the present

paper.

Three versions of NACA 496

The author of the present paper had access to three different versions of NACA 496.  At the conclusion

of their respective main texts, all three versions contain the following lines of text

“LANGLEY MEMORIAL AERONAUTICAL LABORATORY,

     NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS,

          LANGLEY FIELD, VA., May 2, 1934.”
4

indicating that all three are “the same” report.  This common date is important, as will be seen below.

The first version of NACA 496 was in hard‐copy form and was contained in the bound volume

“Twentieth Annual Report of the National Advisory Committee for Aeronautics,” which contains the

significant findings of NACA for fiscal year 1934 (July 1933 through June 1934).  The page numbering

(pages 413 to 433) in this version reflects the fact that NACA 496 was one of many papers in the bound

volume.  Errors in some of the equations, tables, and figures were found (by the present author) to exist

in this version.

The second version was in electronic form and was acquired from the NASA Technical Report Server

website (http://ntrs.nasa.gov).  This version was a reprint of the first version and was issued as a stand‐

alone report.  The date on the cover of this version is 1949.  The page numbering in this version (pages 3

to 23) reflects the fact that it was a stand‐alone report.  Some of the errors found in the first version are

corrected in the second version.

The third version was also a hard‐copy version, was also a stand‐alone report, and was acquired from

the NASA‐Langley technical library.  Pages 3 to 23 in this version are identical to those in the second

version, but the date on the cover of the third version is 1940.

The first tripping point for the present author was the fact that there are three versions of the paper,

each with a different publication date.  A second tripping point was that some of the errors in the first

version were corrected in the second and third versions, but no known errata was issued by NACA, as

confirmed by the staffs of the NASA‐Langley and NASA‐Ames technical libraries.  (Errors from the first

version that were not corrected in the second and third versions, and therefore remain in the second

and third versions, will be discussed throughout the present paper.)

The third version, a corrected version, is contained in Appendix A of the present paper.   

Referring to Equations

NACA 496 contains many equations, but only a minority is identified by letters or numbers.  In the main

text, capital Roman letters and Roman and Arabic numerals, all within parentheses, are used as

identifiers.  In Appendix I of NACA 496, Arabic numerals within parentheses are used as identifiers; in

addition, some terms are identified by Arabic numerals, but without parentheses.  In what would be

very unconventional by today’s report‐writing standards, within NACA 496, there are multiple instances

of different equations being identified using the same equation number.     

In the present paper, references to specific equations in NACA 496 will cite the NACA 496 page and

equation numbers.  For example, equation (A) on page 10, equation (XXI) on page 12, and equation (9)

on page 14 would be identified herein as “equation (10/A),” “equation (12/XXI),” and “equation (14/9),”

respectively.   

Equations in the present paper are identified by Arabic numerals within parentheses.
5

Computational Shortcuts

Given the very rudimentary computational aids of the early 1930’s (paper and pencil, slide rules, four‐

function electro‐mechanical calculators), it is natural to deduce that in the construction of their solution

methods engineers of the time would have sought ways to minimize the overall number and/or to

simplify the complexity of calculations so as to minimize the human time and effort required to obtain a

solution.  It should not be surprising, then, that such time‐ and effort‐saving techniques are evident in

the NACA 496 solution method.  These techniques will be given the general term “computational

shortcuts.”  They will be identified throughout the present paper and are collected in Table III of the

present paper.

III.  OVERVIEW OF NACA REPORT No. 496

NACA 496 lays out the theoretical development of aeroelastic flutter for a typical section with degrees

of freedom in torsion (), aileron deflection (), and vertical deflection (at times referred to in NACA 496

as flexure) (h).  NACA 496 has three appendices:  Appendix I presents a detailed description of the steps

of the implementation of the solution procedure, including tables of numerical values of quantities

required in the solutions; Appendix II presents a number of numerical calculations, including a few

comparisons with experimental data; and Appendix III presents integral evaluations for some of the

velocity potentials.  

Theoretical Development

In NACA 496, the theoretical development begins with four simplifying assumptions:  (1) the flow is

potential and non‐stationary; (2) the “wing” is actually a two‐dimensional typical section with no

thickness and therefore with no airfoil shape; (3) the “wing” motions are sinusoidal and infinitesimal;

and (4) the “wing” has no internal or solid friction, resulting in no internal damping forces.

Next, the circulatory and non‐circulatory velocity potentials are developed.  The magnitude of the

circulation, in the form of Theodorsen’s circulation function, is developed.  After this, the aerodynamic

forces and moments are determined via the chordwise integration of the velocity potentials.  Then, the

aerodynamic forces and moments are combined with the inertia and restraining forces and moments to

produce three second‐order differential equations in the three unknowns  , , and h, equations (10/A),

(10/B), and (10/C), reproduced here:
6

where C(k), Theodorsen’s circulation function, is a complex function of reduced frequency, k.  Equation

(10/A) defines the sum of the moments about the elastic axis; equation (10/B), the sum of the moments

about the aileron hinge; and equation (10/C), the sum of the forces on the entire “wing” in the vertical

direction.  Upon examination of equations (10/A) through (10/C) one sees that k appears only implicitly

gin C(k) while velocity, v, appears multiple times.  

Figure 1 in the present paper is a re‐drawn figure 2 from NACA 496.  It illustrates the definitions and

positive senses of many of the important physical parameters appearing in the equations of motion.

Assumed forms of the unknowns, , , and h. ‐ The assumed forms of the unknowns in equations

(10/A) through (10/C) are introduced in NACA 496 as sine functions of the distance, s, traveled by the

wing from the first vortex element.  Employing the distance formula, s = vt, the sine functions in

complex form may be expressed as:  

௩௜௞݁଴ߙ	ൌ ߙ

௕௧

ߚ ൌ	ߚ଴݁௜ሺ௞௩

௕௧ାఝభሻ

݄ ൌ	݄଴݁௜ሺ௞௩

௕௧ାఝమሻ

where 0, 0, and h0 are the (infinitesimal) amplitudes of , , and h, and ߮ଵ and ߮ଶ are phase angles of

 and h with respect to .  The first and second time derivatives of , , and h are:

ߙሶ ൌ ݅݇ ௩

௕

ߙ  and  ߙሷൌ	െ ቀ݇ ௩

௕

ቁ

ଶ

ߙ

ߚሶ ൌ ݅݇ ௩

௕ ߚ  and  ߚሷൌ	െ ቀ݇ ௩௕

ቁ

ଶ

ߚ

݄ሶ ൌ ݅݇ ௩

௕

݄  and  ݄ሷൌ	െ ቀ݇ ௩

௕

ቁ

ଶ

݄

Substitution of assumed forms into, and normalization of, equations. ‐ Making the substitutions of

equations (1a) through (1c) and their time derivatives into equations (10/A) through (10/C) transforms

the latter equations from three simultaneous differential equations into three simultaneous algebraic

(1a)

(1b)

(1c)
7

equations.  The algebraic equations are then normalized by the quantity  ߢ ቀ

௩

௕

݇ቁ

ଶ

and by their

respective exponentials and amplitudes, resulting in the equations taking the form

ሺܣ௔ఈ ൅	Ωఈܺሻߙ	 ൅ 	ܣ௔ఉߚ ൅ 	ܣ௛௔݄ 		ൌ 		0

ܣ௕ఈߙ 	 ൅ 	ሺܣ௕ఉ ൅	Ω௕ܺሻߚ	 ൅ ܣ௛௕݄ 		ൌ 		0

ܣ௖ఈߙ	൅ 	ܣ௖ఉߚ 	 ൅ 	ሺܣ	൅ ௛௖Ω௛ܺሻ݄		 ൌ 		0

The physical constants of the problem (, a, b, c, x, r, x, r, , , and h) reside in the coefficients,

that is in the A’s and the X terms.   

The quantities Aa, etc. are complex functions of reduced frequency, k, with real parts Ra, etc. and

imaginary parts Ia, etc.  The real parts are defined in equations (12/1) through (12/9) and the imaginary

parts in equations (12/11) through (12/19).  (There is no equation (12/10).  This equation numbering

convention was chosen in NACA 496 so that for any given imaginary part, the second digit in its two‐digit

equation number would be the same as the single‐digit equation number of its corresponding real part.)

The X terms are discussed next.

Products X,X, and hX . ‐  In equations (2a), (2b), and (2c) the products X, X, and hX are

real.  They are derived from quantities   ஼ഀ

ெ௕మ (from eqn. (10/A)),   ஼ഁெ௕మ (from eqn. (10/B)), and  

஼೓

ெ௕ (from

eqn. (10/C)), respectively.  Via the substitution and rearrangement of terms and the use of cancelling

expressions in the numerator and denominator, from page 12 of NACA 496 these products are

ߗఈܺ ൌ	

ఈܥ

ߢଶݒܯଶ݇

ൌ 	 ൬

߱ఈݎఈ

௥ݎ௥߱

൰

ଶ 1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

Ωఉܺ ൌ	

ఉܥ

ߢଶݒܯଶ݇

ൌ 	 ൬

߱ఉݎఉ

௥ݎ௥߱

൰

ଶ 1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

Ω௛ܺ ൌ	

ଶܾ௛ܥ

ߢଶݒܯଶ݇

ൌ 	 ൬ ߱௛

௥ݎ௥߱

൰

ଶ 1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

where, to the right of the second equal sign in each equation, X comprises the two right‐most terms

ܺ ൌ	

1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

and the respective ’s, the remaining terms

Ωఈ ൌ 	 ൬

߱ఈݎఈ

௥ݎ௥߱

൰

ଶ

(3a)

(4)

(5a)

(2a)

(2b)

(2c)

(3c)

(3b)
8

Ωఉ ൌ 	 ൬

߱ఉݎఉ

௥ݎ௥߱

൰

ଶ

Ω௛ ൌ 	 ൬ ߱௛

௥ݎ௥߱

൰

ଶ

The quantities r and rr are a reference frequency and reference length, which may be conveniently

chosen.   

Referring back to the normalization that produced equations (2a) through (2c), a critically important

result of the normalization is that the quantity X has been conveniently isolated from the other terms in

these equations and X is the only quantity in these equations that contains the velocity.  The NACA 496

solution methods take advantage of this result in an ingenious way.  As will be seen immediately below,

X and the ’s are clever artifices created by the author of NACA 496 to enable the determination of the

flutter velocity, vf, and flutter reduced frequency, kf.  At times X and the ’s are treated not as the

known quantities expressed in equations (4) and (5), but rather as parameters.  At other times X and the

’s are treated as the known quantities.   

Solution Methods

The solution of equations (2a), (2b), and (2c) is obtained when their determinant is zero

ቮ

௛௔ܣ ఉ௔ܣ ܺΩఈ ൅	ఈ௔ܣ

௛௕ܣ ܺΩఉ ൅	ఉ௕ܣ ఈ௕ܣ

ܺ௛Ω௖௛ ൅	ܣ ఉ௖ܣ ఈ௖ܣ

ቮ = 0

The form of equation (6), det(A + X) = 0, bears some resemblance to the characteristic polynomial of

the eigenvalue problem, det(A ‐ I) = 0.  The quantity X in eqn. (6) would be analogous to  in the

eigenvalue problem and the quantity , a 3x3 diagonal matrix with diagonal elements , , and h,

would be analogous to the identity matrix.  However, in terms of the solution of equation (6), this is

where its resemblance to the eigenvalue problem ends.

Expanding the determinant in equation (6) yields equation (12/XXI), re‐written here

ΩఈΩఉΩ௛ܺଷ ൅ ൫ΩఈΩఉܣ	൅ ௛௖ΩఉΩ௛ܣ௔ఈ ൅	Ω௛Ωఈܣ௕ఉ൯ܺଶ

൅൫Ωఈܯ௔ఈ ൅	Ωఉܯ௕ఉ ൅	Ω௛ܯ௛௖൯ܺ ൅ ܦ ൌ 0

which is third‐order in X with complex coefficients.  In equation (7), D is the determinant in equation (6)

without the X terms and the M’s are the minors of the determinant without the X terms.  Except for

, , and h, all quantities in equation (7) are functions of reduced frequency, k.   

(5b)

(7)

(6)

(5c)
9

The object of the solution method is to find the values of v and k that cause equation (7) to be satisfied.  

Such values are then identified as vf and kf, the flutter velocity and the flutter reduced frequency,

respectively.

Three‐degree‐of‐freedom solution method. – NACA 496 outlines the method by which the three‐

degree‐of freedom flutter problem (that it names the “general case”) may be solved.   

NACA 496 does not recommend solving directly the third‐order equation with complex coefficients (eqn.

(7)).  Instead, it suggests the computational shortcut of creating two equations, each with real

coefficients, by separating the real and imaginary parts of the original.  The first equation is cubic in X

and is obtained from the real parts of the A’s, M’s, and D; the second is quadratic in X and is obtained

from the imaginary parts.  The right hand side of both equations is zero.  This separation of a complex

equation into its real and imaginary parts eliminates complex arithmetic from the solution method.

At this point in the solution method, the artifice of treating X as a parameter, rather than as a known

quantity, is employed.  The cubic and quadratic equations are each solved for their roots, identified

herein as ܺோభ, ܺோమ, ܺோయ, ܺூభ, and ܺூమ for a large number of reduced frequencies.  The XR’s and XI’s are

then plotted on the same set of axes as functions of the inverse of reduced frequency.  Intersections of

any of the XR’s with any of the XI’s identify the values of X and the values of reduced frequency, which

simultaneously satisfy the cubic and quadratic equations, and therefore also solves the original third‐

order equation with complex coefficients.  These values are the flutter values ‐ Xf and kf.

With Xf and the kf known, the artifice of treating X as a parameter is now abandoned and equation (4) is

employed, which when re‐written to solve for flutter velocity, becomes  

	 ൌ௙ ݒ

1

ඥܺ௙

1

ߢ√

௥ݎ௥ܾ߱

݇௙

where  and b are known constants of the problem and r and rr are chosen reference values.

Although many illustrative examples are presented in Appendix II of NACA 496 none are given for the

general case, almost certainly due to the difficulty in the early 1930’s of solving “by hand” a cubic

equation.

Two‐degree‐of‐freedom solution method. – There are three two‐degree‐of‐freedom subcases available

from the three‐degree‐of‐freedom equations:  torsion – aileron; aileron – flexure; and flexure – torsion.  

NACA 496 lists the subcases in the order just given, but then identifies the first one listed as subcase 3,

the second one listed as subcase 2, and the third one listed as subcase 1.  No reason is offered in NACA

496 for this curious pairing of listing order and subcase number.

NACA 496 presents an ingenious solution method for obtaining the flutter velocity and flutter reduced

frequency for the three two‐degree‐of‐freedom subcases.  This solution method again involves

computational shortcuts and the artifices, but, in yet another curiosity, this solution method differs

fundamentally from the recommended solution method for three degrees of freedom.

(8)
10

Whereas equation (12/XXI) is the characteristic equation obtained from the expansion of the 3x3

characteristic determinant (eqn. (7) in the present paper), equations (12/XXII), (12/XXIII), and (12/XXIV)

are the characteristic equations obtained from the expansion of the three 2x2 minors of the

characteristic determinant.  The reader is referred to Appendix A of the present paper.  

Using the torsion – aileron subcase (subcase 3) as the example, NACA 496 goes through the algebraic

steps of this solution method.  The algebraic steps for subcases 2 and 1 mirror those for subcase 3.  

Equation (12/XXII), a second‐order equation in X with complex coefficients, is the starting point.  NACA

496 again employs the computational shortcut of creating two equations, each with real coefficients, by

separating the real and imaginary parts of the original.  The first equation is quadratic in X and is

obtained from the real parts of Aa, Ab, and Mch; the second is linear in X and is obtained from the

imaginary parts.  The right hand side of both equations is zero.   

At this point the two‐degree‐of‐freedom solution method departs from the three‐degree‐of‐freedom

solution method, but like the latter solution method it also employs the artifices:  it treats both X and 

as parameters, rather than as the known quantities in equations (4) and (5a).  The quadratic and linear

equations and parameters X and  are treated as two equations in two unknowns and are solved

conventionally by substitution.  

From page 13 in NACA 496, the first step is to solve the linear equation for X, yielding

ܺ ൌ	െ

௛௖ܯ

ூ

Ωఈܫ௕ఉ ൅	Ωఉܫ௔ఈ

This expression is then substituted into the quadratic equation in X, thereby eliminating X from that

equation and producing, after several steps, equation (13/XXV), which is quadratic in  and re‐written

here

Ωఈ

ଶ ൫ܯ௛௖

ோ ܫ௕ఉ

ଶ െ ܯ௛௖

ூ ܴ௕ఉܫ௕ఉ൯൅	Ω௔ൣെܯ௛௖ூ ൫ܴ௔ఉܫ௕ఈ ൅ ܫ௔ఉܴ௕ఈ൯	൅ 	2ܯ௛௖ோ ܫ௔ఈܫ௕ఉ൧	

௛௖ܯ	൅

ோ ܫ௔ఈ

ଶ െ	ܯ௛௖

ூ ܴ௔ఈܫ௔ఈ = 0

Equation (10) is then solved for its two roots, Ωఈభand Ωఈమ, for a large number of reduced frequencies.  

With the artifice still in place, but recognizing from equation (5a) that  is defined as the square of a

quantity, only the real positive values of Ωఈభand Ωఈమand their corresponding reduced frequencies are

retained.  A plot of  as a function of 1/k, which contains both Ωఈభand Ωఈమ, is produced.

Next, values of X are found by substituting the values of , just obtained, their corresponding values of

reduced frequency, and the values of MI

ch, Ia, and Ib, at corresponding reduced frequencies, into

equation (9).  NACA 496 chooses to present the quantity X as a function of  rather than as a function

of 1/k.

At this point a new quantity, F, the nondimensional flutter factor, is introduced  

(9)

(10)
11

	ൌ ܨ

1

݇

ඨ1

ܺ

F is also treated as a parameter.  Substituting into equation (11) the X’s and their corresponding reduced

frequencies produces a plot of F as a function of , again containing both Ωఈభand Ωఈమ.

The two curves,  vs. 1/k, and F vs. , represent a family of flutter solutions for a range of ’s.  Each

point on the first curve has a corresponding point on the second and each pair of corresponding points

represents a unique flutter solution.   

To obtain a specific flutter solution from the family of solutions, the artifice is abandoned:  quantities

, X, and F are no longer interpreted as parameters and equations (5a) and (4) are employed.  In

addition, when the expression in equation (4) is substituted for X in equation (11), the flutter factor

becomes

ݒߢ√ 	ൌ ܨ

௥ݎ௥ܾ߱

When solved for flutter velocity, equation (12) becomes

ܨ ൌ௙ ݒ

௥ݎ௥ܾ߱

ߢ√

Thus, once a family of flutter solutions is obtained, the following simple two‐step process for obtaining a

specific flutter solution from the family of solutions is performed:

First, problem‐specific values of , r, r, and rr are substituted into equation (5a), yielding

a problem‐specific value of .  This value of  is located in the plot of  vs. 1/k.  The

value of 1/k corresponding to this value of  then yields, via its inverse, the flutter reduced

frequency, kf.

Second, the problem‐specific value of  is located in the plot of F vs. .  The value of F

corresponding to this value of  and problem‐specific values of b, r, rr, and  are

substituted into equation (13), yielding the flutter velocity, vf.

IV.  CHECKS OF NACA 496 EQUATIONS AND SOLUTION METHOD

Beginning with equations (10/A), (10/B), and (10/C) in the “DIFFERENTIAL EQUATIONS OF MOTION”

section of NACA 496, the present author checked all of the equations in the main text and Appendix I.  

With a few exceptions, discussed below, all equations checked.   

(11)

(12)

(13)
12

Missing 1/b Factors in Expressions on Page 12

In checking the equations in NACA 496, the present author duplicated 12 of the 18 expressions in the

left column on page 12.  However, as they appear in NACA 496, the real expressions in equations (12/3),

(12/6) and (12/9) and the imaginary expressions in equations (12/13), (12/16) and (12/19) each lack a

factor of 1/b compared to what the present author found.  These six expressions make up the entirety

of equation (10/C) and none of the six appears in either equation (10/A) or equation (10/B).   

As it turns out, the missing 1/b factors are of no consequence because the left‐hand side of an equation

whose right‐hand side is zero, such as equation (10/C), may be multiplied by any factor without

changing the equation.  The present author believes that, as a computational shortcut, the author of

NACA 496 multiplied both sides of equation (10/C) by the semi‐chord, b, thus eliminating the 1/b factor

from equations (12/3), (12/6), (12/9), (12/13), (12/16), and (12/19) and thereby relieving engineers of

the time of extra, but unnecessary, multiplications in the hand calculations of these expressions.  

Comparing Real Expressions in Main Text and Appendix I

The real expressions in the left column on page 12 (main text), the real expressions in the right column

on page 14 (Appendix I), and the real expressions in the left column on page 15 (Appendix I), term‐by‐

term, share the same equation numbers.  This assignment of equation numbers is a potential tripping

point because, term‐by‐term, the expressions obviously differ from each other.  As it turns out, the

expressions on page 12 are, term‐by‐term, equal to the sums of the respective expressions on pages 14

and 15.  In checking the equations the present author reproduced all of these expressions.  

Table III on page 17 in Appendix I contains a range of values for the expressions in equations (14/1)

through (14/9).

Comparing Imaginary Expressions in Main Text and Appendix I

The imaginary expressions in the left column on page 12 (main text) correspond to the imaginary

expressions in the right column on page 14 (Appendix I).  A cursory examination of these equations

suggests the identity of corresponding expressions because:  (1) the same symbols for the I’s (Ia, Ib, Ich,

etc.) are used in both places; and (2) the same equation numbers are used in both places.

In checking the equations, the present author reproduced the expressions on page 12.  However, in

comparing the expressions on page 12 with those on page 14, it is seen that each expression on page 12

contains a multiplying factor of 1/k while those on page 14 do not.  As with the omission of the factor of

1/b previously discussed, it turns out that the omission of 1/k is also ultimately of no consequence, and

for the same reason.  As with the 1/b omission, the present author believes that the 1/k omission was

also done deliberately by the author of NACA 496 as another computational shortcut, again relieving

engineers of an extra, but unnecessary, multiplication in the hand calculation of these expressions.   
13

The author of NACA 496 does alert the reader to these differences in the I’s, but via a whisper, not a

shout:  a footnote on the bottom of page 14.  No explanation is given as to why 1/k has been eliminated,

yet the elimination has important consequences in the equations employed in the solution method.  The

differences represented a significant tripping point for the present author, who initially missed the

footnote.   

Table IV on page 17 in Appendix I contains a range of values for the expressions in equations (14/11)

through (14/19) (that is, the expressions without the 1/k factor).

Comparing Equations for Solution Method in Main Text and Appendix I

The NACA 496 solution method is derived beginning at the bottom of the left column on page 12 and

continuing through the end of the main text on page 13.  The procedure for implementing the solution

method is contained in Appendix I on pages 14 through 16.  These parts of NACA 496 will be referred to

as the derivation section and the implementation section, respectively.

The solution method involves solving quadratic equations whose coefficients (a, b, and c) contain

numerous sums and differences of products of the R’s and I’s and other terms.  The other terms are the

real and imaginary parts of the minors of the characteristic determinant (MR

ch, MRa, MRb, and MIch,

MI

a, MIb, and referred to hereinafter as “real minors” and “imaginary minors”), which also contain

sums and differences of products of the R’s and I’s.

As mentioned above, to take advantage of a computational shortcut, all the I’s in the implementation

section had their common 1/k factor eliminated.  For clarity these I’s without the 1/k factor will be called

“new I’s.”  The new I’s and any quantities involving the new I’s (the real and imaginary minors and the

quadratic coefficients a, b, and c) are all affected by the elimination of the 1/k factor.

This situation would have caused no confusion had new symbols (perhaps original symbols with tildes)

been introduced in NACA 496 to represent the new I’s.  The original symbols could have been used in

the derivation section; the new symbols in the implementation section.  However, the author of NACA

496 chose not to do this; he chose to retain the original symbols in the implementation section and in so

doing assigned each symbol two meanings, with each meaning specific to a given section of the paper.

However, with careful attention to section‐specific definitions and by doing a little algebra, it is easy to

prove that the quadratic coefficients in the implementation section are proportional, term‐by‐term, to

the quadratic coefficients in the derivation section, with proportionality constant 1/k2

.  Thus, because

multiplying either equation (both of whose right‐hand sides are zero) by the proper quantity (either k2

or 1/k2

) will produce the other, both equations will yield the same roots, and therefore the same values

of vf and kf.
14

V.  RE‐COMPUTATION OF QUANTITIES EMPLOYED IN NACA 496 SOLUTION METHOD

NACA‐496 contains four tables of quantities needed in the implementation of the NACA‐496 solution

method.  For the engineer of the time, these tables, especially the latter two, would have saved valuable

human computation time.

Constants Resulting from the Integration of Velocity Potentials

Table I in NACA 496 contains constants, Ti, associated with integrations of the velocity potentials.  These

constants are functions of the nondimensional chordwise location of the aileron hinge, c.  At the top of

the right column on page 5 of NACA 496 are the integrals in question and it can be seen that there are

two types of integrals that produce the Ti’s:  those whose limits of integration are c to 1 (integrations

over the aileron only) and those whose limits are ‐1 to 1 (integrations over the entire “wing”).

At the bottom of the right column on the same page are the equations for T1 through T14.  However

Table I lists only 11 of the 14 Ti’s.  The omitted Ti’s (T9, T13, and T14) are also functions of a, the elastic axis

location.

Figure 2 in the present paper contains comparison plots of the Ti’s listed in Table I and present

calculations of the Ti’s.  The open circles correspond to the tabular values, which are shown in NACA 496

for five values of c (‐1, ‐ 1⁄2, 0, 1⁄2, 1).  The solid lines correspond to the present calculations, performed in

Matlab® and computed from c = ‐1 to c = 1 in increments of 0.01.  The plots show there is excellent

agreement between the original results and the present results.

Theodorsen Circulation Function

Table II in NACA 496 contains the real part, F, and the negative of the imaginary part, ‐G, of the

Theodorsen circulation function and their component Bessel functions.  All quantities are functions of

reduced frequency, k, and are shown for 16 values of k.

Figure 3 in the present paper contains comparison plots of F and –G listed in Table II and present

calculations of F and –G.  The horizontal axis is logarithmic to more clearly show (compared to what a

linear axis would show) the comparisons at low values of 1/k.  A logarithmic axis “spreads out” the

values of F and –G at the low values of 1/k and compresses them at the high values of 1/k in such a way

that the distance on the page between 1/k of 0.1 and 0.2 is the same as the distance on the page

between 1/k of 10 and 20. The open circles correspond to the tabular values; the solid lines correspond

to the present results, computed in Matlab® from 1/k = 0.1 to 1/k = 40 in increments of 0.004.   

For values of 1/k of 5 and below there is excellent agreement between the original results and the

present results.  For all values of 1/k except 0.1, the original results and the present results differ by less

than one percent.  For values of 1/k of 10 and above the differences vary from 0.2 to 5.6 percent, with
15

the largest differences occurring at 10 and 40.  These largest differences will produce corresponding

differences seen in the illustrative examples.

Reduced‐Frequency‐Dependent Portions of the Elements of the Equations of Motion

Table III in NACA 496 contains the frequency‐dependent portions of the real parts of the coefficients of

the equations of motion; Table IV contains the imaginary parts.  All quantities in Tables III and IV are

functions of reduced frequency, k.  Some quantities in Tables III and IV are functions of the

nondimensional chordwise location of the elastic axis, a; other quantities are functions of the

nondimensional chordwise location of the aileron hinge, c; some quantities are functions of both; and

two quantities are functions of neither.  The quantities in Tables III and IV will be discussed together.  

Figure 4 in the present paper contains comparison plots of the real and imaginary parts from Tables III

and IV and from present calculations, again performed in Matlab®.  As was the case for figure 3, in this

figure the horizontal axis is logarithmic to more clearly show (compared to what a linear axis would

show) the comparisons at low values of 1/k.  Figure 4(a) contains elements from the pitching moment

equation; figure 4(b), from the hinge moment equation; and figure 4(c), from the vertical force

equation.  For the present calculations the inverse of reduced frequency ranged from 1/k = 0.1 to 1/k =

40, in increments of 0.0004.   

The plots are for values of a and c of ‐0.4 and 0.5, respectively.  Except for the comparison plot for ܴ௔௛

ᇱᇱ ,

which will be discussed in the next paragraphs, all comparisons show reasonably good agreement.  The

comparisons at values of 1/k of 10, 20, and 40 consistently show differences and are due to the

differences in F and –G, discussed above.

The comparison plot for ܴ௔௛	

ᇱᇱ in figure 4(a) shows significant differences between the original and

present results.  Close examination of the values in Table III reveals that, for a = ‐0.4, ܴ௔௛	

ᇱᇱ is very nearly

directly proportional to 1/k, that is, is very nearly linear with 1/k.  To more easily see this linearity,

ܴ௔௛	

ᇱᇱ as a function of 1/k has been redrawn in figure 4(d) with a linear horizontal axis.  A red dashed line

has been added to emphasize the linear relationship between ܴ௔௛	

ᇱᇱ and 1/k.

If one assumes the values of ܴ௔௛	

ᇱᇱ in Table III are correct and then employs equation (14/3) to back

calculate to find the values of G necessary to produce these values of ܴ௔௛

ᇱᇱ , one finds that the resulting

G’s are independent of 1/k and equal to the value of G at 1/k = 1.  (According to Table II in NACA 496

and figure 3 in the present paper, G varies from about ‐0.2 to 0 and clearly is a function of 1/k.)  The ܴ௔௛

ᇱᇱ

row of Table III corresponding to a = ‐0.4 and c = 0.5 is in error.  

VI.  RE‐COMPUTATION OF ILLUSTRATIVE EXAMPLES IN NACA 496

This section of the paper presents the re‐computation of the illustrative examples contained in Appendix

II of NACA 496.  These examples comprise all three two‐degree‐of‐freedom subcases that may be
16

derived from the three‐degree‐of‐freedom general case.  NACA 496 presents no examples for the

general case.

A brief orientation to figures 5 through 17 in the present paper:  figure‐by‐figure they contain scanned

reproductions of the original figures 5 through 17 from NACA 496, but with the addition of the present

re‐computations.  In all cases, the original results are in black and the re‐computed results are in color.  

Because of “fuzziness” resulting from the scans, the figure legends for these figures have been re‐typed.

Additionally, in figures 5 through 11 and 15 through 17, the colored  and  symbols are results from

the present implementation, but using the tabular values in Tables I, II, III, and IV of NACA 496, and the

colored solid and dashed lines are results from the present implementation using present computed

values throughout.   

As mentioned in Section III of the present paper, the plots of  vs 1/k and F vs  are comprised of two

values of , the two quadratic roots obtained from the solution of equation (13/XXV) (as well as two

values each from the solutions of eqns. (13/XXVI) and (13/XXVII)).  With this in mind, the colored solid

lines and the colored  symbols correspond to the first quadratic root; the colored dashed lines and the

colored  symbols correspond to the second quadratic root.

Implementation of NACA 496 Solution Method

The implementation of the NACA 496 solution method was performed by the present author via the

writing and execution of m‐files in Matlab®.  The equations outlined in the main text and Appendix I of

NACA 496 were solved for 100,001 values of the inverse of reduced frequency, beginning at zero, with

an increment of 0.0005.  

Standard Case

Figures 5 through 10 present what NACA 496 calls a “standard” case, represented by the quantities:

 = 0.1;  c = 0.5;  a = ‐0.4;  x = 0.2;  r

2 = 0.25;  x = 1⁄80;  r2 = 1⁄160;  and , , h variable.

Figures 5 and 6 present results for Subcase 3, torsion‐aileron; figures 7 and 8 present results for Subcase

2, aileron‐flexure; and figures 9 and 10 present results for Subcase 1, flexure‐torsion.   

For these three subcases, the first figure in each pair shows its respective plot of  versus 1/k; the

second figure in each pair shows its respective plot of F versus .  However, the quantity F shown in

these figures is not the quantity F defined on the bottom of page 15 and the top of page 16 of NACA 496

(and by eqn. (11) in the present paper).  The quantity F shown in figures 6, 8, and 10 is actually the

square of the defined quantity.  This double use of a symbol to represent not only its original definition,

but also the square of its original definition, represented a significant tripping point for the present

author.
17

Subcase 3. – Figures 5 and 6 contain the results for this subcase.  While there are some differences

between the original and present results in these figures, the overall magnitudes and shapes of the

respective curves are very similar.  In figure 5, the agreement between original and present results is

excellent at values of 1/k greater than 0.5, but progressively less so at lower values of 1/k.  In figure 6,

the agreement is better at the lower values of  than at the higher values.  In both figures, the present

results represented by the colored symbols and colored lines exhibit excellent agreement with each

other, reflecting the excellent agreement between the original tabular quantities and the present

computed quantities, as illustrated in figure 4.

Subcase 2. ‐ Figures 7 and 8 contain the results for this subcase.  In each figure, there are three sets of

results, labeled (a), (b), and (c), representing three different values of x, the center of gravity of the

aileron.  The original and present curves (a) and (b) in figure 7 agree fairly well in the overall magnitudes

and shapes of the respective curves.  However, the original and present curves (c) do not agree well in

either magnitude or shape.  The original and present curves in figure 8 echo the agreement and

disagreement noted for the curves in figure 7.  In both figures, the present results represented by the

colored symbols and colored lines exhibit excellent agreement with each other, reflecting again the

excellent agreement between the original tabular quantities and the present computed quantities, as

illustrated in figure 4.

Subcase 1. ‐ Figures 9 and 10 contain the results for this subcase.  Once again, in these figures the

original and present results are very similar in terms of the overall magnitudes and shapes of the

respective curves.  In figure 9, the original curve does not “touch” the vertical axis, whereas the present

curve does (as did all original and present curves in figures 5 and 7).  This trend is repeated in figure 10,

but for the horizontal axis.  The present results represented by the colored symbols and colored lines do

not agree as well in figure 9 and 10 as they did in figures 5 and 6 and figures 7 and 8.   These differences

are due in part to the errors in ܴ௔௛

ᇱᇱ in Table III, and illustrated in figure 4, as previously discussed.

Parameter Variations for Subcase 1

Figures 11 through 14 contain plots of flutter factor, F, as defined on pages 15 and 16 of NACA 496, as

functions of various parameters.  Figure 11 presents results for variation in frequency ratio, 1/2;

figure 12, for variation in the location of the elastic axis, a; figure 13, for variation in the radius of

gyration, r; and figure 14, for variation in the location of the center of gravity, x.  All results in this

subsection of the paper are for Subcase 1, flexure‐torsion.   

For figures 12, 13, and 14, as a check on the present results, independent calculations were performed

by a colleague of the present author who implemented in Matlab® the equations found in reference 5

for a typical section, including unsteady aerodynamics defined by the Theodorsen circulation function.   

The code loops on velocity.  For each velocity, the calculations are run through an iterative loop until the

frequencies are converged, followed by an evaluation of damping values to determine if the system is

stable or unstable.   
18

To be able to directly compare the independent flutter results with the present flutter results, the

independent flutter velocities were converted to flutter factors using equation (12).  As will be seen

below, in all cases, the present results and the independent results are in excellent agreement, differing

in F by no more than 0.5% and in reduced frequency by no more than 3.7%, thereby simultaneously

confirming and affording confidence in the present results.  

Variation of Frequency Ratio. – The original figure 11 contains theoretical and experimental results, in

black.  The present results are shown in colored symbols and colored curves.  The quantities specific to

this example are:

 = 1⁄400;  a = ‐0.4;  x = 0.2;  r = 0.5;  1/2 variable.

This figure is similar to figure 10, but as already stated, the F in this figure is the as‐defined F.  In

addition, in this figure, equation 5(c) has been employed to convert the units of the abscissa from h to

1/2 (that is, h/).

Across the limited range of frequency ratio chosen for this figure, the original theoretical results and

both sets of present results are in excellent agreement.  The colored symbols in the figure correspond to

values of 1/k of 10, 20, and 40.  The re‐computed curve in this figure does not “touch” the horizontal

axis as did the re‐computed curve in figure 10 because the much smaller value of mass ratio, , for

figure 11 changes the character of the solution.

Variation of Elastic Axis Location. – In addition to containing the dependency of F on elastic axis

location, a, figure 12 also contains the dependency of a quantity D on a.  According to the sign

convention established in figure 2 of NACA 496, increasingly negative values of a corresponds to a

forward movement of the elastic axis.

The quantities specific to this example are:

 = 1⁄4;  1/2 = 1⁄6;  x = 0.2;  r = 0.5;  a variable.  

In the present effort to re‐compute D as a function of a, another tripping point was encountered:  D is

never defined in NACA 496.  However, the present author attempted to deduce the definition of D so

that a comparison could be made between the original and the re‐computed curves.   

In the left column on page 16 of NACA 496, there is a short discussion that ends with the equations (un‐

numbered in NACA 496 but numbered herein)

ோݒ	ൌ஽ ݒ

ଵ

భ

మ

ା௔    

ܨோݒ	ൌ ிݒ

where a and F are defined as above, vF is the flutter velocity, vR is a reference velocity defined to be the

divergence velocity for the specific case of the elastic axis located at the mid‐chord (that occurs when a

= 0), and vD is the divergence velocity for other locations of the elastic axis.   

(14)

(15)
19

Velocity vR is a proportionality constant in equations (14) and (15).  In equation (15), vR relates the flutter

velocity, vF, to the flutter factor, F.  The present author speculated that in equation (14) vR might relate

the divergence velocity, vD, to a “divergence factor,” D.  If this were so, then according to equation (14),

quantity D would be equal to the inverse of (1⁄2 + a).  As a test of this hypothesis the inverse of (1⁄2 + a)

was plotted as a function of a.  The hypothesis failed; the curve did not match curve D in figure 12.   

(Equation (14) is clearly incorrect because, according to the definition of vR, when a is zero, vD should be

equal to vR, but under this condition equation (14) predicts vD to be twice vR.)

This miss‐match led the present author to re‐derive the expression relating vD to D, which resulted not in

equation (14), but in the following expression  

ோݒ	ൌ஽ ݒ

1

√1 ൅ 2ܽ

When the quantity ଵ

√ଵାଶ௔

is plotted as a function of a, the colored dashed curve in figure 12 results.  

Except for an area of slight disagreement (less than 5%) between a = ‐0.2 and a = ‐0.4, the original and

present curves agree with each other very well.  This overall agreement led the present author to

conclude that the original curve D in figure 12 represents the expression the author of NACA 496 knew

to be correct (eqn. (16)) and that the expression on page 16 of NACA 496 (equivalent to eqn. (14)) is a

misprint.

The colored circles in figure 12 represent the re‐calculation of curve F.  Each circle was obtained by

selecting a value of a and executing the NACA 496 solution method for Subcase 1 in which a curve

similar to that in figure 10 is used to find the value of F at the following value of h ‐ 

ߗ ௛ൌ 	 ቀ ఠ೓

ఠഀ௥ഀ

ቁ

ଶ

=  ቀ ఠభ

ఠమ௥ഀ

ቁ

ଶ

=  ቀ

ఠభ

ఠమ

ቁ

ଶ

ቀ ଵ

௥ഀ

ቁ

ଶ

ൌ	 ቀ

ଵ

଺

ቁ

ଶ

ቀ ଵ

଴.ହቁ

ଶ

=  1⁄9

It is obvious that the re‐calculation differs markedly from the original.  Not only do the numerical values

significantly differ from each other, but the trends with increasingly negative values of a are opposite

each other.  The three large X’s at values of a of ‐0.2, ‐0.3, and ‐0.4 are from the independent calculation

and are in excellent agreement with the present calculation.  The present author can offer no

explanation for the significant disagreement between the present calculation and the original.

Variation of Radius of Gyration. – The original figure 13 illustrates the dependency of F on the radius of

gyration, r. for two different values of frequency ratio, 1/2.  Curve A corresponds to 1/2 = 1⁄6; curve

B corresponds to 1/2 = 1.  Other quantities specific to this example are

 = 1⁄4;  a = ‐0.4;  x = 0.2;  r variable

For ease of presentation and in an attempt to preserve the figure numbering convention described at

the beginning of this section, this figure is divided into “offspring” figures:  13A and 13B.  Further, figure

13B has parts (a), (b), and (c).

(16)
20

Figure 13A contains colored circles representing the re‐calculations intended to match the original curve

A.  Each circle was obtained by selecting a value of r and executing the NACA 496 solution method for

Subcase 1.  The agreement between the original and present results is very good at values of r of 0.75

and greater, but slightly less so at lower values of r.  If original curve A is extrapolated to the left, the

difference between the original and present values of F at r = 0.25 is about 7.5%.  The large X’s at

values of r of 0.50, 0.75, and 1.00 are from the independent calculation and are in excellent agreement

with the present calculations.

Figure 13B(a) contains colored squares representing the re‐calculations intended to match the original

curve B.  Each square was also obtained by selecting a value of r and executing the NACA 496 solution

method for Subcase 1.  The trends of the original and present results agree (F decreases with increasing

values of r), but in terms of absolute numbers the agreement is very poor:  at the higher values of r,

the magnitudes of the original results are more than a factor of two larger than those of the present

results; at the lower values of r they approach a factor of two smaller.  In fact, the vertical scale of the

original figure had to be extended in order to accommodate the present result at r = 0.25.   

N.B.:  In the following discussion there are two uses of the word “square.”  The first use refers to the

product of a quantity and itself; the second use refers to the geometric shape of the colored plotting

symbols in figure 13B.   

In contemplating the poor agreement between the original and present results, the present author

noticed that the values of the present results, represented by the colored squares, at r = 0.375 and 1.0

(about 1.5 and 0.25, respectively) are approximately equal to the squares of the values of the original

curve at the same r’s (about 1.25 and 0.5, respectively).  This relationship led the present author to

compute and plot the squares of the values of the original curve B (curve labeled in the figure as the

“square of curve B”).   

As can be seen in figure 13B(b), the agreement between the present results and the square of curve B is

much better than the agreement between the present results and curve B itself.  The better agreement

may be quantified by examining the sums of the squared differences between the results:  those of the

former are less than half those of the latter.

This better agreement suggests that an error occurred in the creation of the original curve B.  Because

the author of NACA 496 uses the symbol F to mean, at times, a quantity and, at other times, the square

of that quantity, there is ample room for one meaning to be used when the other meaning was

intended.  The present author believes such an error occurred in the creation of curve B (but not in the

creation of curve A).

In figure 13B(c), the large X’s at values of r of 0.50, 0.75, and 1.00 are from the independent calculation

and are in excellent agreement with the present calculations and is supporting evidence that the original

curve B was incorrectly plotted.
21

Variation of Center of Gravity Location. – The original figure 14 illustrates the dependency of F on the

location of the center of gravity, x, for three different combinations of quantities.  Each combination is

represented by a different curve, labeled A, B, and C.  According to the sign convention established in

figure 2 of NACA 496, increasingly positive values of x corresponds to an aft movement of the center of

gravity.

Because of the complexity of the results to be presented for this figure and in an attempt to preserve

the figure numbering convention described at the beginning of this section, figure 14 will be divided up

into “offspring” figures:  14A, 14B, and 14C.  In addition, each offspring figure will have parts (a), (b), and

(c).   

According to the figure legend in the original figure 14, the quantities specific to this example are

for curve A:  = 1⁄400;  a = ‐0.4;  1/2 = 1⁄6;  r = 0.5;  x variable

and

for curve B:  = 1⁄4;  a = ‐0.4;  1/2 = 1⁄6;  r = 0.5;  x variable

and

for curve C:  = 1⁄100;  a = ‐0.4;  1/2 = 1;  r = 0.5;  x variable.

The colored circles in Figures 14A are the result of present calculations intended to match the original

curve A; the colored squares in Figures 14B, to match the original curve B; and the colored diamonds in

Figures 14C, to match the original curve C.  Each present result represented by a symbol was obtained

by selecting a value of x and executing the NACA 496 solution method for Subcase 1.   

The original figure 14 contains conflicting information regarding curve C.  The figure legend defines the

quantities listed above; but, within the figure itself is a notation indicating other quantities – “1/2 = 1,  

4 = .01” – with a leader pointing to curve C.  The ratio of frequencies is the same as the ratio specified

in the figure legend; however, the value of  (= 1⁄400) differs from that in the figure legend.  Present

calculations were performed for both values of .

In figure 14A(a), the original results and the present results differ significantly, showing opposite trends

with increasing values of x.  However, the present author observed the following:  the present value of

F at x = 0 is very close to the original value of F at x = 0.4 and vice versa; the present value of F at x =

0.1 is very close to the original value of F at x = 0.3 and vice versa; and the present value of F at x = 0.2

is very close to the original value of F at x = 0.2.

This observation led to figure 14A(b), in which the original curve A has been rotated about the vertical

axis x = 0.2.  The rotated curve A and the colored circles exhibit excellent agreement.

In Figure 14A(c) the large X’s at values of x of 0.1, 0.2, and 0.3 are from the independent calculation and

are in excellent agreement with the present results and the rotated curve A.  This triple agreement

suggested to the present author that original curve A was incorrectly plotted.   
22

In the early 1930’s, intricate and multi‐stepped engineering calculations (such as those required to

produce the results in NACA 496) were performed by hand; results were recorded and transcribed by

hand; and then plotted by hand.  It is not unrealistic to consider that with the many steps performed by

humans, errors occurred.  The present author believes that in the case of curve A, transcription errors

occurred – correct values of F were paired with non‐corresponding values of x – which then led to an

incorrect plot.

The trends observed in figures 14A are repeated in figures 14B and 14C:  in figures 14B(a) and 14C(a) the

original results and the present results (represented now by the colored squares and diamonds) differ

significantly, showing opposite trends with increasing values of x; in figures 14B(b) and 14C(b) the

original curves B and C have been rotated about the vertical axis x = 0.2 producing excellent agreement

with the present results; and in figures 14B(c) and 14C(c) the large X’s from the independent calculation

are in excellent agreement with the present results and the rotated curves B and C.  Again, this time for

curves B and C, the present author believes that transcription errors occurred – correct values of F were

paired with non‐corresponding values of x – which then led to incorrect plots.

Wings A and B

The final three figures in NACA 496 (figs. 15, 16, and 17) present results for two wings, designated in

NACA 496 as Wing A and Wing B.  For each wing, experimental results were obtained.  Both wings had

symmetrical airfoils with chords of 12.7 centimeters (5 inches) and thickness‐to‐chord ratios of 0.12.  

Both wings were tested at zero degrees incidence.   

The text of Appendix II of NACA 496 states that Wing A was constructed of aluminum with quantities:  

 = 1∕416;  a = ‐0.4;  x = 0.31, 0.173, and 0.038;  r

2 = 0.33;   = 7x2.

The text of Appendix II of NACA 496 states that Wing B was constructed of wood, with a trailing‐edge

control surface, with quantities:

 = 1⁄100;  c = 0.5;  a = ‐0.4;  x = 0.192;  r

2 = 0.178;   

x = 0.019;  r

2 = 0.0079;   = 17.6x2.

Figure 15 contains Subcase‐1 results for Wing A in which flutter velocity is plotted as a function of

frequency ratio h /.

Figure 16 contains Subcase‐2 results for Wing B in which flutter velocity is plotted as a function of

frequency ratio  /h.

Figure 17 contains Subcase‐3 results for Wing B in which flutter velocity is plotted as a function of

frequency ratio /.
23

To obtain the curves in these figures, first the appropriate subcase‐specific NACA 496 solution method

was executed, resulting in curves of subcase‐specific F (as defined) vs. subcase‐specific .  Next, for each

value of F equation (13) was employed, converting the ordinates in these curves from nondimensional F

to dimensional vf.  Finally, the equations for the subcase‐specific  ’s from equations (5a), (5b), and (5c)

were employed, converting the abscissas from subcase‐specific   to subcase‐specific frequency ratio.

Referring back to equation (13), it can be seen that to obtain the dimensional flutter velocity reference

values r and rr are required, which differ from subcase to subcase.  For Subcase 1, values of  and r

are required; for Subcase 2, a value of h is required; and for Subcase 3, values of  and r  are

required.  NACA 496 supplies values for , r, and r, but, in another tripping point, it does not supply

values for h and .  Therefore, for the present calculations for Subcases 2 and 3, values of h and 

had to be chosen.  The choice was based on a trial‐and‐error investigation of candidate values and the

selection criterion for the chosen values of h and  was that there be reasonable agreement between

the original and the present curves.   

Subcase 1. – Figure 15 contains theoretical and experimental results, in black.  The present results are

shown in colored  and  symbols and solid and dashed colored curves.  Even though three values of x

are listed above in the quantities for this subcase, x = 0.173 was chosen for the present calculation

because it produced the best agreement with the original curve in figure 15.  Over the entire frequency

range the trends for original theoretical results and both sets of present results are in very good

agreement.  The overall agreement over the first half of the range is excellent with maximum

differences less than 5 percent; over the second half, less so with maximum differences of about 15

percent.

Subcase 2. – Figure 16 contains theoretical and experimental results, in black.  The present results are

shown in colored  and  symbols and solid and dashed colored curves.  Theoretical results are given

for two values of quantity xb:  x = 0.019 and x = 0.01.  For this subcase, an assumed value of h of

5.8x2 was selected because it produced reasonable agreement between the original results and the

present results.   

In examining the colored curves, it is seen that the upper (corresponding to higher values of v) and

lower (lower values of v) extremes of both curves extend to a frequency ratio of zero.  Both colored

curves are comprised of solid and dashed portions, indicating that two roots of equation (10) are

present in each solution.  The transition from solid to dashed occurs between frequency ratios of 0.8

and 0.9.  The calculations that produced these curves had a range of 1/k from 0 to 500.

In examining the black theoretical curves, it is seen that the upper extremes of both curves terminate at

a frequency ratio of about 0.75 and the lower extremes terminate at a frequency ratio of about 0.82.  

The difference in horizontal extent between the black and colored curves suggests to the present author

that the calculations that produced the original curves had a very limited range of 1/k.  Further, because

the frequency ratio at the lower extremes of the black curves corresponds roughly to the frequency

ratio at the transition of the colored curves from solid to dashed, it is possible that the black curves only

capture a single root of equation (13/XXVI).
24

Subcase 3. – Figure 17 contains theoretical and experimental results, in black.  The present results are

shown in colored  and  symbols and solid and dashed colored curves.  For this subcase, an assumed

value of  of 4.4x2 was selected because it produced reasonable agreement between the original

results and the present results.   

As was observed in figure 16 for subcase 2 because the frequency ratio at the lower extreme of the

black curve corresponds roughly to the frequency ratio at the transition of the colored curve from solid

to dashed, it is possible that the black curve only captures a single root of equation (13/XXV).

VII.  CONCLUDING REMARKS

This paper has reported on an extensive study of reference 2, NACA Report No. 496, “General Theory of

Aerodynamic Instability and the Mechanism of Flutter,” by Theodore Theodorsen (referred to herein as

NACA 496), including the re‐computation of its numerical results.  The study included reading and re‐

reading the paper, checking its equations, and checking its solution methods.  With the few exceptions

noted in the present paper, the study revealed the equations and solution methods of NACA 496 to be

correct.  Because of the crude (by today’s standards) computational tools available at the time NACA

496 was written, the author of NACA 496 conceived of and employed several time‐ and effort‐saving

computational shortcuts, pointed out and summarized in a table in the present paper.  The study also

revealed some “tripping points,” aspects of NACA 496 that have the potential to cause confusion, and

identified some errors in NACA 496 that were not corrected in the later versions of the paper.  Tripping

points and errors are also summarized in tables in the present paper.    

Appendix II in NACA 496 contains a number of numerical examples, all of which were re‐computed in

the present paper.  The re‐computed results were overlaid on the original figures from NACA 496.  The

re‐computations were accomplished by employing the method developed in the original work, but by

using modern computational tools.  In most instances the overall shapes of the original and re‐

computed curves are similar, with varying degrees of agreement in magnitude.  In some instances, there

is significant dissimilarity between the original and re‐computed curves.  Most dissimilarities are

attributed to human computational and transcription errors that are speculated to have occurred during

the creation of NACA 496.  One dissimilarity (namely, that regarding the quantity F in fig. 12) is

unattributed.  Several of the re‐computed results were checked by independent calculations and, in

each instance, the re‐computed and independent results were in excellent agreement.  

Appendix A of the present paper contains NACA 496.

Appendix B contains a Matlab® program that implements the two‐degree‐of‐freedom flutter equations

and the two‐degree‐of‐freedom solution method found in Appendix I of NACA 496.  By executing this

program, results in figures 5 through 11 and 15 through 17, found in Appendix II of NACA 496, may be

re‐computed directly.  This program produces both physical and non‐physical solutions.  Physical

solutions are based on real values of the quantity F and the ’s; non‐physical solutions are based on

complex values.
25

Appendix C of the present paper offers three alternate solution methods, all of which derive from the

two‐degree‐of‐freedom flutter equations of NACA 496.  Compared to the number of computations

involved in the NACA 496 solution method, these alternate solution methods are much more

computationally intensive and would have been excluded from consideration in 1934 for that reason.  

To accuracies within a tenth of a percent, all three alternate solution methods yielded the same value of

flutter velocity, vf, and the same value of flutter reduced frequency, kf, as those predicted by the

solution method of NACA 496.   

Appendix D of the present paper offers a solution of the three‐degree‐of‐freedom flutter equations

following the method proposed in NACA 496.  This solution method produced an example result for the

“standard case.”  In addition, for systematic variations in the dimensional modal frequencies, in the

limit, the three‐degree‐of‐freedom results (vf and kf) approached the two‐degree‐of‐freedom results,

including the number of flutter modes predicted.
26

APPENDIX A

NACA REPORT NO. 496
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

APPENDIX B

MATLAB® PROGRAM FOR RE‐COMPUTING RESULTS IN NACA 496

This appendix contains a Matlab® program that implements the two‐degree‐of‐freedom flutter

equations and the two‐degree‐of‐freedom solution method found in Appendix I of NACA 496.  By

executing this program, results in figures 5 through 11 and 15 through 17, found in Appendix II of NACA

496, may be re‐computed directly.  The re‐computation of results in figures 12 through 14 is discussed

below.

The program is structured as depicted in the simplified flow chart, below, on the left.  There are four

“groups” of re‐computed figures that may be chosen, according to the table, below, on the right.  The

user inputs an integer from 1 to 4, thereby specifying a group.  Then, based on the integer, the program

selects a set of quantities (, x, r, x, r, etc.), computes the coefficients of the equations of motion,

solves the equations of motion and plots the results.  The “cases” in the table are consistent with

subsection titles in section VI of the main body of the present paper.

  

The re‐computation of results in figures 12 through 14 requires multiple executions of the program, with

the user selecting group number 2 and supplying appropriately modified sets of quantities (see figure

legends of figs. 12 through 14) for this group.

During the normal execution of the solution method, the computed quantity F and the computed ’s

may be, at times, real and, at other times, complex.  Only real values of F and the ’s represent physical

solutions; complex values of F and the ’s represent non‐physical solutions.   

The Matlab® program in this appendix produces both physical and non‐physical solutions.  When

plotting, Matlab® does not distinguish between real quantities and complex quantities and, if a quantity

is complex, Matlab® will plot its real part and issue a warning message in the command window.  The

user is cautioned to check the command window for such messages, and if present, carefully examine

the variable F and the ’s to discern which values in the solution are real and which are complex.  By so
53

doing the user may determine if a particular locus of points in a figure represents a physical solution or a

non‐physical solution.     

‐‐‐‐‐‐‐‐  ‐‐‐‐‐‐‐‐  ‐‐‐‐‐‐‐‐ 

%% This program reproduces figs 5 --> 11 and 15 --> 17 from NACA 496

% Figs 12 --> 14 require multiple executions of group number 2 with

% appropriately modified input quantities

%

% Implemented using NACA 496 equations and 2DOF solution method

% Boyd Perry, III

% NASA-Langley Research Center

% May 2015

clc

clear

close all

%% Define k and kinv

N = 10001;

oneoverkmax = 100;

deltaoneoverk = oneoverkmax/(N-1);

for i=1:N

 kinv(i) = (i-1)*deltaoneoverk;

 k(i) = 1/kinv(i);

end

%% Define Quantities

a = [0 -0.2 -0.4];

c = [0 0.5];

% For all calculations a = -0.4 and c = 0.5

group = input('Enter group number (1, 2, 3, or 4) ')

% Group number may be 1, 2, 3, or 4

% Group number 1 is for the Standard Case and reproduces figs 5 --> 10

% Group number 2 is for the Standard Case and reproduces fig 11

% Group number 3 is for Wing A and reproduces fig 15

% Group number 4 is for Wing B and reproduces figs 16 & 17

% Input quantities for group number 1

if group == 1

kappa = 0.1;

xa = 0.2;

ra2 = 0.25;

xb = 1/80; % Corresponds to curve (a) in figs 7 & 8

%xb = 1/40; % Corresponds to curve (b) in figs 7 & 8

%xb = 1/160; % Corresponds to curve (c) in figs 7 & 8

rb2 = 1/160;

end
54

% Input quantitie for group number 2

if group == 2

kappa = 1/400;

xa = 0.2;

ra2 = 0.25;

xb = 1/80;

rb2 = 1/160;

end

% Input quantitie for group number 3

if group == 3

kappa = 1/416;

%xa = 0.31;

xa = 0.173;

%xa = 0.038;

ra2 = 0.33;

xb = 1/80;

rb2 = 1/160;

omegaa = 7*2*pi;

b = 2.5/12;

end

% Input quantitie for group number 4

if group == 4

kappa = 1/100;

xa = 0.192;

ra2 = 0.178;

xb = 0.019;

%xb = 0.01;

rb2 = 0.0079;

omegaa = 17.6*2*pi;

omegah = 5.8*2*pi;

omegab = 4.4*2*pi;

b = 2.5/12;

end

% Other relationships

ra = sqrt(ra2);

rb = sqrt(rb2);

%% Define array indices

%

% For Ti constants:

% Ti(kk)

% kk ~ variation of quantity c

%

% For R' arrays:

% R'(j,kk)

% j ~ variation of quantity a

% kk ~ variation of quantity c

%

% For F and G, real and imag parts of Theodorsen Function:

% F(i) and G(i)
55

% i ~ variation of reduced frequency

%

% For R" and I arrays:

% R"(i,j,kk) and I(i,j,kk)

% i ~ variation of reduced frequency

% j ~ variation of quantity a

% kk ~ variation of quantity c

%

%% Define Ti

for kk=1:2 % loop on quantity c [0, 0.5]

 c2(kk) = c(kk)*c(kk);

 TiTerm1 = 1 + c2(kk);

 TiTerm2 = 1 - c2(kk);

 TiTerm3 = sqrt(TiTerm2);

 TiTerm4 = acos(c(kk));

 TiTerm5 = 1/8 + c2(kk);

 TiTerm6 = 7 + 2*c2(kk);

 T1(kk) = -(1/3)*TiTerm3*(2+c2(kk)) + c(kk)*TiTerm4;

 T2(kk) = c(kk)*TiTerm2 - TiTerm3*TiTerm1*TiTerm4 + c(kk)*TiTerm4*TiTerm4;

 T3(kk) = -TiTerm5*TiTerm4*TiTerm4 + 

(1/4)*c(kk)*TiTerm3*acos(c(kk))*TiTerm6 - (1/8)*TiTerm2*(5*c2(kk)+4);

 T4(kk) = -TiTerm4 + c(kk)*TiTerm3;

 T5(kk) = -TiTerm2-TiTerm4*TiTerm4 + 2*c(kk)*TiTerm3*TiTerm4;

 T6(kk) = T2(kk);

 T7(kk) = -TiTerm5*TiTerm4 + (1/8)*c(kk)*TiTerm3*TiTerm6;

 T8(kk) = -(1/3)*TiTerm3*(2*c2(kk)+1) + c(kk)*TiTerm4;

 T10(kk) = TiTerm3 + TiTerm4;

 T11(kk) = TiTerm4*(1-2*c(kk)) + TiTerm3*(2-c(kk));

 T12(kk) = TiTerm3*(2+c(kk)) - acos(c(kk))*(2*c(kk)+1);

end

%% Define F and G

for i=1:N

 J0(i) = besselj(0,k(i));

 J1(i) = besselj(1,k(i));

 Y0(i) = bessely(0,k(i));

 Y1(i) = bessely(1,k(i));

 TheoTerm1(i) = J1(i) + Y0(i);

 TheoTerm2(i) = Y1(i) - J0(i);

 

 Fnumer(i) = J1(i)*TheoTerm1(i) + Y1(i)*TheoTerm2(i);

 Gnumer(i) = Y1(i)*Y0(i) + J1(i)*J0(i);

 

 denom(i) = TheoTerm1(i)^2 + TheoTerm2(i)^2;

 

 F(i) = Fnumer(i)/denom(i);

 G(i) = -Gnumer(i)/denom(i);

end

%% Define R' Terms

for kk=1:2 % loop on quantity c [0, 0.5]

 for j=1:3 % loop on quantity a [0 -0.2 -0.4]
56

Raap(j,kk) = -ra2/kappa - (1/8 + a(j)*a(j));

Rabp(j,kk) = -rb2/kappa - (c(kk)-a(j))*xb/kappa + T7(kk)/pi + (c(kk)-

a(j))*T1(kk)/pi;

Rahp(j,kk) = -xa/kappa + a(j);

Rbap(j,kk) = Rabp(j,kk);

Rbbp(j,kk) = -rb2/kappa + T3(kk)/(pi*pi);

Rbhp(j,kk) = -xb/kappa + T1(kk)/pi;

Rcap(j,kk) = Rahp(j,kk);

Rcbp(j,kk) = Rbhp(j,kk);

Rchp(j,kk) = -1/kappa - 1;

 end

end

%% Define R" and I Terms

% Equations (1)->(9) and (11)->(19) on p. 14

for kk=1:2 % loop on quantity c [0, 0.5]

 p = -(1/3)*(1-c(kk)*c(kk))^(3/2);

 for j=1:3 % loop on quantity a [0 -0.2 -0.4]

 for i=1:N % loop on reduced frequency

 

 RITerm1 = (a(j)+0.5);

 RITerm2 = (0.5-a(j))*G(i) - kinv(i)*F(i);

 RITerm3 = T11(kk)*G(i)-2*kinv(i)*T10(kk)*F(i);

 

 Raapp(i,j,kk) = kinv(i)*2*RITerm1*RITerm2;

 Rabpp(i,j,kk) = 

kinv(i)*(1/pi)*((T4(kk)+T10(kk))*kinv(i)+RITerm1*RITerm3);

 Rahpp(i,j,kk) = kinv(i)*2*RITerm1*G(i);

 

 Rbapp(i,j,kk) = -kinv(i)*(T12(kk)/pi)*RITerm2;

 Rbbpp(i,j,kk) = -kinv(i)*(1/pi)*(1/pi)*((T12(kk)/2)*RITerm3-

kinv(i)*(T5(kk)-T4(kk)*T10(kk)));

 Rbhpp(i,j,kk) = -kinv(i)*(T12(kk)/pi)*G(i);

 

 Rcapp(i,j,kk) = -kinv(i)*2*RITerm2;

 Rcbpp(i,j,kk) = -kinv(i)*(1/pi)*RITerm3;

 Rchpp(i,j,kk) = -kinv(i)*2*G(i);

 

 RITerm4 = (0.5-a(j))*F(i) + kinv(i)*G(i);

 RITerm5 = T11(kk)*F(i)+2*kinv(i)*T10(kk)*G(i);

 

 Iaa(i,j,kk) = -2*RITerm1*RITerm4 + (0.5-a(j));

 Iab(i,j,kk) = -(1/pi)*(RITerm1*RITerm5 + 2*p + (0.5-

a(j))*T4(kk));

 Iah(i,j,kk) = -2*RITerm1*F(i);

 Iba(i,j,kk) = (T12(kk)/pi)*RITerm4 +(1/pi)*(p-T1(kk)-0.5*T4(kk));

 Ibb(i,j,kk) = (0.5/(pi*pi))*(T12(kk)*RITerm5-T4(kk)*T11(kk));

 Ibh(i,j,kk) = (T12(kk)/pi)*F(i);
57

 Ica(i,j,kk) = 2*RITerm4 + 1;

 Icb(i,j,kk) = (1/pi)*(RITerm5-T4(kk));

 Ich(i,j,kk) = 2*F(i);

 

 end

 end

end 

%% Create R's and I's for j=3 (a=-0.4) and kk=2 (c=0.5)

j = 3;

kk = 2;

for i=1:N

 

 R1(i) = Raap(j,kk) + Raapp(i,j,kk);

 R2(i) = Rabp(j,kk) + Rabpp(i,j,kk);

 R3(i) = Rahp(j,kk) + Rahpp(i,j,kk);

 

 R4(i) = Rbap(j,kk) + Rbapp(i,j,kk);

 R5(i) = Rbbp(j,kk) + Rbbpp(i,j,kk);

 R6(i) = Rbhp(j,kk) + Rbhpp(i,j,kk);

 

 R7(i) = Rcap(j,kk) + Rcapp(i,j,kk);

 R8(i) = Rcbp(j,kk) + Rcbpp(i,j,kk);

 R9(i) = Rchp(j,kk) + Rchpp(i,j,kk);

 

 I11(i) = Iaa(i,j,kk);

 I12(i) = Iab(i,j,kk);

 I13(i) = Iah(i,j,kk);

 I14(i) = Iba(i,j,kk);

 I15(i) = Ibb(i,j,kk);

 I16(i) = Ibh(i,j,kk);

 I17(i) = Ica(i,j,kk);

 I18(i) = Icb(i,j,kk);

 I19(i) = Ich(i,j,kk);

end

%% Solve Case 3 Example Using Equations from Appendix I (Figs 5 & 6)

if group == 1

 

% Compute quadratic coefficients

for i=1:N

 

MRch(i) = R1(i)*R5(i) - R2(i)*R4(i) - kinv(i)*kinv(i)*(I11(i)*I15(i)-

I12(i)*I14(i));

MIch(i) = R1(i)*I15(i) - R2(i)*I14(i) + I11(i)*R5(i) - I12(i)*R4(i);
58

A(i) = MRch(i)*I15(i)*I15(i) - MIch(i)*R5(i)*I15(i);

B(i) = -MIch(i)*(R2(i)*I14(i)+ I12(i)*R4(i)) + 2*MRch(i)*I11(i)*I15(i);

C(i) = MRch(i)*I11(i)*I11(i) - MIch(i)*R1(i)*I11(i);

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);

 OmegaAlpha1(i) = r(1);

 OmegaAlpha2(i) = r(2);

 X1(i) = -MIch(i)/(OmegaAlpha1(i)*I15(i)+I11(i));

 X2(i) = -MIch(i)/(OmegaAlpha2(i)*I15(i)+I11(i));

 

 F1(i) = kinv(i)*kinv(i)/X1(i); % This is actually the square of F

 F2(i) = kinv(i)*kinv(i)/X2(i); % This is actually the square of F

end

figure(5)

plot(kinv, OmegaAlpha1, 'r.', kinv, OmegaAlpha2, 'b.');

grid

axis([0 1 0 160])

xlabel('1/k','FontWeight','bold');

ylabel('Omega-Alpha','FontWeight','bold');

title('Fig 5, Case 3, Standard Case','FontWeight','bold');

figure(6)

plot(OmegaAlpha1, F1, 'r.', OmegaAlpha2, F2, 'b.');

grid

axis([0 180 0 22])

xlabel('Omega-Alpha','FontWeight','bold');

ylabel('F','FontWeight','bold');

title('Fig 6 Case 3, Standard Case','FontWeight','bold');

end

%% Solve Case 2 Example Using Equations from Appendix I (Figs 7 & 8)

if group == 1

% Compute quadratic coefficients

for i=1:N

 

MRaa(i) = R5(i)*R9(i) - R6(i)*R8(i) - kinv(i)*kinv(i)*(I15(i)*I19(i)-

I16(i)*I18(i));

MIaa(i) = R5(i)*I19(i) - R6(i)*I18(i) + I15(i)*R9(i) - I16(i)*R8(i);

A(i) = MRaa(i)*I19(i)*I19(i) - MIaa(i)*R9(i)*I19(i);

B(i) = -MIaa(i)*(R6(i)*I18(i)+ I16(i)*R8(i)) + 2*MRaa(i)*I15(i)*I19(i);

C(i) = MRaa(i)*I15(i)*I15(i) - MIaa(i)*R5(i)*I15(i);
59

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);

 OmegaBeta1(i) = r(1);

 OmegaBeta2(i) = r(2);

 X1(i) = -MIaa(i)/(OmegaBeta1(i)*I19(i)+I15(i));

 X2(i) = -MIaa(i)/(OmegaBeta2(i)*I19(i)+I15(i));

 

 F1(i) = kinv(i)*kinv(i)/X1(i); % This is actually the square of F

 F2(i) = kinv(i)*kinv(i)/X2(i); % This is actually the square of F

end

figure(7)

plot(kinv, OmegaBeta1, 'r.', kinv, OmegaBeta2, 'b.');

grid

axis([0 4 -0.008 0.014])

xlabel('1/k','FontWeight','bold');

ylabel('Omega-Beta','FontWeight','bold');

title('Fig 7, Case 2, Standard Case','FontWeight','bold');

annotation('textbox','String',{['xb = ', num2str(xb)]},'FontWeight','bold', 

'FontSize',10,...

 'BackgroundColor',[1 1 1]);

figure(8)

plot(OmegaBeta1, F1, 'r.', OmegaBeta2, F2, 'b.');

grid

axis([0 0.014 0 1.5])

xlabel('Omega-Beta','FontWeight','bold');

ylabel('F','FontWeight','bold');

title('Fig 8, Case 2, Standard Case','FontWeight','bold');

annotation('textbox','String',{['xb = ', num2str(xb)]},'FontWeight','bold', 

'FontSize',10,...

 'BackgroundColor',[1 1 1]);

end

%% Solve Case 1 Example Using Equations from Appendix I (Figs 9 & 10)

if group == 1

% Compute quadratic coefficients

for i=1:N

 

MRbb(i) = R9(i)*R1(i) - R7(i)*R3(i) - kinv(i)*kinv(i)*(I19(i)*I11(i)-

I17(i)*I13(i));

MIbb(i) = R9(i)*I11(i) - R7(i)*I13(i) + I19(i)*R1(i) - I17(i)*R3(i);

A(i) = MRbb(i)*I11(i)*I11(i) - MIbb(i)*R1(i)*I11(i);
60

B(i) = -MIbb(i)*(R7(i)*I13(i)+ I17(i)*R3(i)) + 2*MRbb(i)*I19(i)*I11(i);

C(i) = MRbb(i)*I19(i)*I19(i) - MIbb(i)*R9(i)*I19(i);

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);

 OmegaH1(i) = r(1);

 OmegaH2(i) = r(2);

 X1(i) = -MIbb(i)/(OmegaH1(i)*I11(i)+I19(i));

 X2(i) = -MIbb(i)/(OmegaH2(i)*I11(i)+I19(i));

 

 F1(i) = kinv(i)*kinv(i)/X1(i); % This is actually the square of F

 F2(i) = kinv(i)*kinv(i)/X2(i); % This is actually the square of F

end

figure(9)

plot(kinv, OmegaH1, 'r.', kinv, OmegaH2, 'b.');

grid

axis([0 10 -30 120])

xlabel('1/k','FontWeight','bold');

ylabel('Omega-H','FontWeight','bold');

title('Fig 9, Case 1, Standard Case','FontWeight','bold');

figure(10)

plot(OmegaH1, F1, 'r.', OmegaH2, F2, 'b.');

grid

axis([0 22 0 1.5])

xlabel('Omega-H','FontWeight','bold');

ylabel('F','FontWeight','bold');

title('Fig 10, Case 1, Standard Case','FontWeight','bold');

end

%% Solve Case 1 Example Using Equations from Appendix I (Fig 11)

if group == 2

 

% Compute quadratic coefficients

for i=1:N

 

MRbb(i) = R9(i)*R1(i) - R7(i)*R3(i) - kinv(i)*kinv(i)*(I19(i)*I11(i)-

I17(i)*I13(i));

MIbb(i) = R9(i)*I11(i) - R7(i)*I13(i) + I19(i)*R1(i) - I17(i)*R3(i);

A(i) = MRbb(i)*I11(i)*I11(i) - MIbb(i)*R1(i)*I11(i);

B(i) = -MIbb(i)*(R7(i)*I13(i)+ I17(i)*R3(i)) + 2*MRbb(i)*I19(i)*I11(i);

C(i) = MRbb(i)*I19(i)*I19(i) - MIbb(i)*R9(i)*I19(i);
61

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);

 OmegaH1(i) = r(1); 

 OmegaH2(i) = r(2); 

 X1(i) = -MIbb(i)/(OmegaH1(i)*I11(i)+I19(i)); 

 X2(i) = -MIbb(i)/(OmegaH2(i)*I11(i)+I19(i)); 

 

 F1(i) = kinv(i)/sqrt(X1(i)); % F as defined on pp. 14-15

 F2(i) = kinv(i)/sqrt(X2(i)); % F as defined on pp. 14-15

 

 om1om21(i) = ra*sqrt(OmegaH1(i)); 

 om1om22(i) = ra*sqrt(OmegaH2(i)); 

end

figure(11)

plot(om1om21, F1, 'r.', om1om22, F2, 'b.');

grid

axis([0 1.6666666667 0 1.75])

xlabel('omega-1 / omega-2','FontWeight','bold');

ylabel('F','FontWeight','bold');

title('Fig 11, Case 1, Standard Case (but with kappa = 

1/400)','FontWeight','bold');

end

%% Solve Case 1 Example Using Equations from Appendix I (Fig 15)

if group == 3

 

% Compute quadratic coefficients

for i=1:N

 

MRbb(i) = R9(i)*R1(i) - R7(i)*R3(i) - kinv(i)*kinv(i)*(I19(i)*I11(i)-

I17(i)*I13(i));

MIbb(i) = R9(i)*I11(i) - R7(i)*I13(i) + I19(i)*R1(i) - I17(i)*R3(i);

A(i) = MRbb(i)*I11(i)*I11(i) - MIbb(i)*R1(i)*I11(i);

B(i) = -MIbb(i)*(R7(i)*I13(i)+ I17(i)*R3(i)) + 2*MRbb(i)*I19(i)*I11(i);

C(i) = MRbb(i)*I19(i)*I19(i) - MIbb(i)*R9(i)*I19(i);

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);
62

 OmegaH1(i) = r(1);

 OmegaH2(i) = r(2);

 X1(i) = -MIbb(i)/(OmegaH1(i)*I11(i)+I19(i));

 X2(i) = -MIbb(i)/(OmegaH2(i)*I11(i)+I19(i));

 

 F1(i) = kinv(i)/sqrt(X1(i)); % F as defined on pp. 14-15

 F2(i) = kinv(i)/sqrt(X2(i)); % F as defined on pp. 14-15

 

 fac1 = b*omegaa*ra/sqrt(kappa);

 fac2 = 12/39.37; % convert from fps to mps

 

 vel1(i) = F1(i)*fac1*fac2;

 vel2(i) = F2(i)*fac1*fac2;

 

 om1om21(i) = ra*sqrt(OmegaH1(i));

 om1om22(i) = ra*sqrt(OmegaH2(i));

end

figure(15)

plot(om1om21, vel1, 'r.', om1om22, vel2, 'b.');

grid

axis([0 1.5 0 50])

xlabel('omega-h / omega-alpha','FontWeight','bold');

ylabel('Velocity, meters per sec','FontWeight','bold');

title('Fig 15, Case 1, Wing A','FontWeight','bold');

annotation('textbox','String',{['xa = ', num2str(xa)]},'FontWeight','bold', 

'FontSize',10,...

 'BackgroundColor',[1 1 1]);

end

%% Solve Case 2 Example Using Equations from Appendix I (Fig 16)

if group == 4

 

% Compute quadratic coefficients

for i=1:N

 

MRaa(i) = R5(i)*R9(i) - R6(i)*R8(i) - kinv(i)*kinv(i)*(I15(i)*I19(i)-

I16(i)*I18(i));

MIaa(i) = R5(i)*I19(i) - R6(i)*I18(i) + I15(i)*R9(i) - I16(i)*R8(i);

A(i) = MRaa(i)*I19(i)*I19(i) - MIaa(i)*R9(i)*I19(i);

B(i) = -MIaa(i)*(R6(i)*I18(i)+ I16(i)*R8(i)) + 2*MRaa(i)*I15(i)*I19(i);

C(i) = MRaa(i)*I15(i)*I15(i) - MIaa(i)*R5(i)*I15(i);

end

% Solve quadratic equation

for i=2:N

 Quad = [A(i) B(i) C(i)];
63

 r = roots(Quad);

 OmegaBeta1(i) = r(1);

 OmegaBeta2(i) = r(2);

 X1(i) = -MIaa(i)/(OmegaBeta1(i)*I19(i)+I15(i));

 X2(i) = -MIaa(i)/(OmegaBeta2(i)*I19(i)+I15(i));

 

 F1(i) = kinv(i)/sqrt(X1(i)); % F as defined on pp. 14-15

 F2(i) = kinv(i)/sqrt(X2(i)); % F as defined on pp. 14-15

 fac1 = b*omegah/sqrt(kappa);

 fac2 = 12/39.37; % convert from fps to mps

 

 vel1(i) = F1(i)*fac1*fac2;

 vel2(i) = F2(i)*fac1*fac2;

 

 om1om21(i) = sqrt(OmegaBeta1(i))/rb;

 om1om22(i) = sqrt(OmegaBeta2(i))/rb;

end

figure(16)

plot(om1om21, vel1, 'r.', om1om22, vel2, 'b.');

grid

axis([0 1.8 0 50])

xlabel('omega-beta / omega-h','FontWeight','bold');

ylabel('Velocity, meters per sec','FontWeight','bold');

title('Fig 16, Case 2, Wing B','FontWeight','bold');

annotation('textbox','String',{['xb = ', num2str(xb)]},'FontWeight','bold', 

'FontSize',10,...

 'BackgroundColor',[1 1 1]);

end

%% Solve Case 3 Example Using Equations from Appendix I (Fig 17)

if group == 4

 

% Compute quadratic coefficients

for i=1:N

 

MRch(i) = R1(i)*R5(i) - R2(i)*R4(i) - kinv(i)*kinv(i)*(I11(i)*I15(i)-

I12(i)*I14(i));

MIch(i) = R1(i)*I15(i) - R2(i)*I14(i) + I11(i)*R5(i) - I12(i)*R4(i);

A(i) = MRch(i)*I15(i)*I15(i) - MIch(i)*R5(i)*I15(i);

B(i) = -MIch(i)*(R2(i)*I14(i)+ I12(i)*R4(i)) + 2*MRch(i)*I11(i)*I15(i);

C(i) = MRch(i)*I11(i)*I11(i) - MIch(i)*R1(i)*I11(i);

end

% Solve quadratic equation

for i=2:N
64

 Quad = [A(i) B(i) C(i)];

 r = roots(Quad);

 OmegaAlpha1(i) = r(1);

 OmegaAlpha2(i) = r(2);

 X1(i) = -MIch(i)/(OmegaAlpha1(i)*I15(i)+I11(i));

 X2(i) = -MIch(i)/(OmegaAlpha2(i)*I15(i)+I11(i));

 

 F1(i) = kinv(i)/sqrt(X1(i)); % F as defined on pp. 14-15

 F2(i) = kinv(i)/sqrt(X2(i)); % F as defined on pp. 14-15

 fac1 = b*omegab*rb/sqrt(kappa);

 fac2 = 12/39.37; % convert from fps to mps

 

 vel1(i) = F1(i)*fac1*fac2;

 vel2(i) = F2(i)*fac1*fac2;

 

 om1om21(i) = sqrt(OmegaAlpha1(i))*(rb/ra);

 om1om22(i) = sqrt(OmegaAlpha2(i))*(rb/ra);

end

figure(17)

plot(om1om21, vel1, 'r.', om1om22, vel2, 'b.');

grid

axis([0 2.8 0 40])

xlabel('omega-alpha / omega-beta','FontWeight','bold');

ylabel('Velocity, meters per sec','FontWeight','bold');

title('Fig 17, Case 3, Wing B','FontWeight','bold');

end

 
65

APPENDIX C

ALTERNATE TWO‐DEGREE‐OF‐FREEDOM SOLUTION METHODS

Three alternate solution methods are offered, all of which are based on the two‐degree‐of‐freedom

flutter equations of NACA 496.  Compared to the number of computations involved in the NACA 496

solution method, these alternate solution methods are much more computationally intensive and would

have been excluded from consideration in 1934 for that reason.  But with today’s computational

resources, these alternate solution methods require only a few seconds of execution time on a desktop

computer.  The solution methods will be presented in order of increasing complexity.

Although these alternate solution methods are appropriate for all three subcases mentioned in NACA

496, the equations for Subcase 1, flexure – torsion, are chosen to illustrate the methods.  The starting

point for each alternate solution method is equation (XXIV) on page 12 of NACA 496, quadratic in X with

complex coefficients, re‐written here   

Ω௛Ωఈܺଶ ൅ ሺΩ௛ܣ௔ఈ ൅ Ωఈܣ௛௖ሻܺ ൅	ܯ௕ఉ ൌ 0

with companion equations

ܺ ൌ	

1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

Ω௛ ൌ 	 ቀ ఠ೓

ఠೝ௥ೝ

ቁ

ଶ

Ωఈ ൌ ൬

߱ఈݎఈ

௥ݎ௥߱

൰

ଶ

where, for Subcase 1, the reference quantities, r and rr, are chosen to be  and r, respectively.

An illustrative example will be presented for each alternate solution method.  The problem‐specific

quantities chosen for the example are those employed in NACA 496 for its “standard case”:

 = 0.1;  c = 0.5;  a = ‐0.4;  x = 0.2;  r

2 = 0.25;  

and with  = 100, h = 50, and b = 1

Using the appropriate values in equations (C3) and (C4), the values of h and  are both 1.   

Because equation (C1) is quadratic, zero, one, or two flutter solutions are possible.  However, for the

problem‐specific quantities chosen, only one flutter solution is obtained.

As expected and as will be seen, all three alternate solution methods will yield almost exactly the same

value of flutter velocity, vf, and almost exactly the same value of flutter reduced frequency, kf, as those

predicted by the two‐degree‐of‐freedom solution method of NACA 496.  The family of solutions for

Subcase 1 for the standard case already exists in figures 9 and 10.  When h = 1 is used in these figures

(C3)

(C4)

(C1)

(C2)
66

the problem‐specific solution is obtained:  vf = 173.26 feet per second (fps) and kf = 0.4355.  The first

row of Table CI contains this result.

On occasion throughout this appendix, the abbreviation ASM will be used to denote “alternate solution

method.”

Alternate Solution Method No. 1

This alternate solution method is directly analogous to the three‐degree‐of‐freedom solution method

discussed in NACA 496.  It employs the computational shortcut of creating two equations, each with real

coefficients, by separating the real and imaginary parts of equation (C1), thereby eliminating complex

arithmetic.  The right‐hand side of both new equations is zero.   

Ω௛Ωఈܺଶ ൅ ሺΩ௛ܴ௔ఈ ൅ Ωఈܴ௖௛ሻܺ ൅	ܯ௕ఉ

ோ ൌ 0

ሺΩ௛ܫ௔ఈ ൅ Ωఈܫ௛௖ሻܺ ൅	ܯ௕ఉ

ூ ൌ 0

The first equation, (C5), is quadratic in X and is obtained from the real parts of Aa, Ach, and Mb; the

second equation, (C6), is linear in X and is obtained from the imaginary parts.   

At this point in ASM No. 1, the artifice of treating X as a parameter, rather than as a known quantity, is

employed.  Equations (C5) and (C6) are each solved for their roots, identified herein as ܺோభ, ܺோమ and ܺூ,

for a large number of reduced frequencies.  But, because from equation (C2), it is seen that X is

proportional to the inverse square of velocity, only the real positive values of ܺோభ, ܺோమ and ܺூ are

retained.   

Next, the real positive roots are plotted as functions of the inverse of reduced frequency on the same

set of axes.  The intersection of either of the XR’s with the XI identifies the value of X and the value of

reduced frequency that simultaneously solve equations (C5) and (C6), and therefore also solves the

original quadratic equation with complex coefficients, equation (C1).  These values are the flutter values

‐ Xf and kf.  

With Xf and kf known, the artifice of treating X as a parameter is abandoned and equation (C2) is

employed, which when re‐written to solve for flutter velocity, becomes  

	 ൌ௙ ݒ

ଵ

ඥ௑೑

ଵ

√఑

௕ఠഀ௥ഀ

௞೑

Plugging Xf, kf, and the other quantities into equation (C7) yields the flutter velocity.

Illustrative Example for Alternate Solution Method No. 1. – ASM No. 1 was implemented in Matlab®

where equations (C5) and (C6) were solved 10,000 times for values of the inverse of reduced frequency

ranging from 0.001 to 10 in increments of 0.001.  Figure C1 presents the results from this

implementation: the ordinate represents the real positive roots of equations (C5) and (C6); the abscissa

(C7)

(C5)

(C6)
67

is the inverse of reduced frequency.  The solid black line represents the locus of the ܺோభ, the first real

positive root of equation (C5); the dashed black line represents the locus of the ܺோమ, the second real

positive root of equation (C5); and the solid red line represents the locus of the ܺூ, the real positive root

of equation (C6).  The circle represents the intersection of the ܺோభ locus with the ܺூ locus and was

obtained using an interpolation scheme.   

By inverting the horizontal coordinate of the intersection, the flutter reduced frequency is obtained.  By

using it and the vertical coordinate in equation (C7), the flutter velocity is obtained.  Both results are

contained in the second row of Table CI and are seen to agree extremely well with the NACA 496 results.

Alternate Solution Method No. 2

This alternate solution method is closely related to ASM No. 1 and also eliminates complex arithmetic

via the same computational shortcut.  But, rather than initially treating X as a parameter, as was done in

ASM No. 1, in this solution method, X is assigned values and the flutter solution is obtained via

systematic trial and error.   

ASM No. 2 employes an outer loop on velocity and an inner loop on reduced frequency and, using

equation (C2), it computes a value of X at each passing.  These X’s are then substituted into the

following equations  

Ω௛Ωఈܺଶ ൅ ሺΩ௛ܴ௔ఈ ൅ Ωఈܴ௖௛ሻܺ ൅	ܯ௕ఉ

ோ ൌ ܼଵ

ሺΩ௛ܫ௔ఈ ൅ Ωఈܫ௛௖ሻܺ ൅	ܯ௕ఉ

ூ ൌ ܼଶ

producing values of Z1 and Z2 at each passing.   

Recall that the flutter condition is obtained when Z1 and Z2 are both zero for the same values of v and k.  

However, unless luck intervenes or infinitesimally small increments in v and k are chosen for the outer

and inner loops, it is almost assured that no combination of v and k will be found, which results in Z1 and

Z2 being exactly zero simultaneously.  This being the case, two different interpolation schemes are

employed to find the combination of v and k that does result in Z1 and Z2 being exactly zero

simultaneously.

When each inner loop is complete (that is, for each velocity), curves of Z1 and Z2 as functions of reduced

frequency are created.  Because equation (C8) is quadratic, there can be zero, one, or two values of

reduced frequency that cause Z1 to be zero.  Because equation (C9) is linear, there can be zero or one

value of reduced frequency that causes Z2 to be zero.   

The Z1 and Z2 curves may or may not pass through zero.  For curves that do pass through zero, a first

interpolation scheme is employed to identify the values of reduced frequency at the zero crossings and

these values and their corresponding velocity are potential solutions and are stored in arrays.  There can

(C8)

(C9)
68

be up to three arrays:  two for the Z1’s and one for the Z2’s.  These identified reduced frequencies are

defined as ݇௭భభ

, ݇௭భమ and ݇௭మ, respectively.  

When the outer loop is complete, the reduced frequencies at the zero crossings, ݇௭భభ

, ݇௭భమand ݇௭మ,  are

plotted as functions of velocity.  Each point on the ݇௭భభ

and 	݇௭భమcurves (that is, each reduced frequency

and its companion velocity) corresponds to Z1 being exactly equal to zero.  Likewise, each point on the

݇௭మ curve corresponds to Z2 being exactly equal to zero.  The flutter condition is obtained when either

the ݇௭భభ

curve or the	݇௭భమcurve intersects the ݇௭మcurve.  Via a second interpolation scheme, the values of

v and k at the intersection are determined.  These values are vf and kf because they cause Z1 and Z2 to be

zero simultaneously.

Illustrative Example for Alternate Solution Method No. 2. – This alternate solution method was also

implemented in Matlab® where equations (C2), (C8), and (C9) were solved within nested loops.  The

inner loop varied reduced frequency from 0.01 to 2 in increments of 0.01; the outer loop varied velocity

from 160 fps to 185 fps in increments of 0.5 fps.  (The limits of the outer loop were chosen based on the

knowledge of flutter velocity gained from the previous illustrative example.)

Figure C2 presents a plot of Z1 and Z2 as functions of reduced frequency for a velocity of 160 fps, well

below the flutter velocity.  As can be seen  Z1, the solid curve, intersects zero at two points, ݇௭భభ

= 0.3706

and ݇௭భమ

= 0.5317, and Z2, the dashed curve, intersects zero at one point, ݇௭మ= 0.4740.  These values of

reduced frequency are combined with similar reduced frequencies for the other velocities to form the

arrays that are plotted in the next figure.

Figure C3 presents a plot of ݇௭భభ

, ݇௭భమ,	and , ݇௭మas functions of velocity.  As can be seen ݇௭భమ, the dashed

curve, intersects the ݇௭మ, the dashed‐dotted curve, at the values of velocity and reduced frequency

indicated in the third row of Table CI.  As with the results for ASM No. 1, those for ASM No. 2 are in

excellent agreement with the results from NACA 496.

Alternate Solution Method No. 3

This alternate solution method is very straightforward, but uses complex arithmetic.  It also initially

employs the artifice of treating X as a parameter and later abandons it.   

Treating X as a parameter, equation (C1), quadratic in X with complex coefficients, is solved for X directly

for a large number of reduced frequencies.  Because the coefficients are complex the resulting X’s are

also complex, but the X’s are not complex conjugates.  These complex X’s are designated X1 and X2.

With X1 and X2 and their corresponding reduced frequencies known, the artifice is abandoned and

equation (C10)

	ൌ ݒ

1

√ܺ

1

ߢ√

ܾ߱ఈݎఈ

݇

(C10)
69

is employed, yielding velocities, v1 correspondng to X1, and v2, corresponding to X2, that are also

complex and also not complex conjugates.     

All the complex X’s, and therefore all the complex v’s and their corresponding k’s, satisfy equation (C1)

and are therfore all flutter solutions.  However, because in actual practice the flutter velocity is real, not

complex, the task now is to examine all of the complex v’s to determine if any have a zero imaginary

part, that is to determine if any are real.   

A convenient way to perform this examination is to plot the imaginary part of v1 against the real part of

v1 and the imaginary part of v2 against the real part of v2, yielding two curves in the complex plain.  Each

point on each curve corresponds to a different value of reduced frequency.  If one of the curves crosses

the real axis, at the crossing, the imaginary part is zero, the real part is the flutter velocity, and the

examination is complete.  Interpolation is then employed, taking the complex velocities immediately

before and immediately after the crossing to determine the value of the corresponding real part at the

crossing, which is vf, and the corresponding value of reduced frequency at the crossing, which is kf.  

Illustrative Example for Alternate Solution Method No. 3. – This alternate solution method was also

implemented in Matlab® where equations (C1) and (C10) were solved 10,000 times for reduced

frequencies ranging from 0.005 to 50 in increments of 0.005.  

Figure C4 presents plots of the imaginary parts of v1 and v2 as functions of their respective real parts.  

The solid line is v1; the dashed line, v2.  By the nature of the curves, it is apparent that, as stated above,

v1 and v2 are not complex conjugates.  As can be seen in the figure, the dashed curve crosses the real

axis and the interpolated values of velocity and reduced frequency at the crossing are indicated in the

fourth row of Table CI.  The results for ASM No. 3 are in perfect agreement with the results from NACA

496.

Table CI

Comparison of Flutter Predictions for

Standard Case

    

Solution Method vf kf

   fps    

    

NACA 496 173.26 0.4355

    

Alt. Soln. Meth. No. 1 173.22 0.4356

    

Alt. Soln. Meth. No. 2 173.21 0.4358

    

Alt. Soln. Meth. No. 3 173.26 0.4355

   
70

APPENDIX D

THREE‐DEGREE‐OF‐FREEDOM SOLUTION METHOD

This appendix provides a solution to the three‐degree‐of‐freedom flutter equation following the

procedure recommended in NACA 496.  Equation (XXI) on page 12 of NACA 496, re‐written here as

equation (D1), is the three‐degree‐of‐freedom flutter equation

ΩఈΩఉΩ௛ܺଷ ൅ ൫ΩఈΩఉܣ	൅ ௛௖ΩఉΩ௛ܣ௔ఈ ൅	Ω௛Ωఈܣ௕ఉ൯ܺଶ

൅൫Ωఈܯ௔ఈ ൅	Ωఉܯ௕ఉ ൅	Ω௛ܯ௛௖൯ܺ ൅ ܦ ൌ 0

with companion equations

ܺ ൌ	

1

൬ ߢ

௥߱௥ݎܾ

൰ ݇ݒ

ଶ

Ωఈ ൌ ቀ

ఠഀ௥ഀ

ఠೝ௥ೝ

ቁ

ଶ

Ωఉ ൌ ൬

߱ఉݎఉ

௥ݎ௥߱

൰

ଶ

Ω௛ ൌ 	 ቀ ఠ೓

ఠೝ௥ೝ

ቁ

ଶ

In equation (D1) the A’s, M’s, and D are complex quantities, making equation (D1) cubic in X with

complex coefficients.

The subsection entitled “Solution Methods” within section III of the main body of the present paper

outlines the three‐degree‐of‐freedom solution method.   To recap, this solution method involves the

artifice of treating X as a parameter and the computational shortcut of creating two equations, each

with real coefficients, by separating the real and imaginary parts of equation (D1), shown here as

equations (D6) and (D7)

ΩఈΩఉΩ௛ܺଷ ൅ ൫ΩఈΩఉܴ௖௛ ൅	ΩఉΩ௛ܴ௔ఈ ൅	Ω௛Ωఈܴ௕ఉ൯ܺଶ

൅൫Ωఈܯ௔ఈ

ோ ൅	Ωఉܯ௕ఉ

ோ ൅	Ω௛ܯ௛௖

ோ ൯ܺ ൅ ܦோ ൌ 0

൫ΩఈΩఉܫ	൅ ௛௖ΩఉΩ௛ܫ௔ఈ ൅	Ω௛Ωఈܫ௕ఉ൯ܺଶ

൅൫Ωఈܯ௔ఈ

ூ ൅	Ωఉܯ௕ఉ

ூ ൅	Ω௛ܯ௛௖

ூ ൯ܺ ൅ ܦூ ൌ 0

Equations (D6) and (D7) are each solved for their roots, identified herein as ܺோభ, ܺோమ, ܺோయ, ܺூభ and ܺூమ, for

a large number of reduced frequencies.  But, because from equation (D2), it is seen that X is

proportional to the inverse square of velocity, only the real positive values of ܺோభ, ܺோమ, ܺோయ, ܺூభ and

ܺூమare of interest.   

(D1)

(D5)

(D6)

(D7)

(D2)

(D3)

(D4)
71

Next, all of the real positive roots are plotted as functions of the inverse of reduced frequency on the

same set of axes.  The intersection(s) of any of the XR’s with any of the XI’s identify(ies) the value(s) of X

and the value(s) of reduced frequency that simultaneously solve equations (D6) and (D7), and therefore,

also solves the original cubic equation with complex coefficients, equation (D1).  These values are the

flutter values ‐ Xf and kf.

With Xf and the kf known, the artifice of treating X as a parameter is abandoned and equation (D2) is

employed, which when re‐written to solve for flutter velocity, becomes  

	 ൌ௙ ݒ

1

ඥܺ௙

1

ߢ√

௥ݎ௥ܾ߱

݇௙

Plugging Xf, kf, and the other quantities into equation (D8) yields the flutter velocity.

The three‐degree‐of‐freedom solution method is capable of finding zero, one, two, or three flutter

modes, depending on the problem‐specific quantities chosen.

Illustrative Example for Three‐Degree‐of‐Freedom Solution Method

The problem‐specific quantities chosen for this example are those employed in NACA 496 for its

“standard case”:

 = 0.1;  c = 0.5;  a = ‐0.4;  x = 0.2;  r

2 = 0.25;  x = 1⁄80;  r2 = 1⁄160;

and with  = 100;   = 125;  h = 50;  b = 1;  r = 1;  rr = 1

This solution method was implemented in Matlab® where equations (D6) and (D7) were solved 10,000

times for values of the inverse of reduced frequency ranging from 0.0005 to 5 in increments of 0.0005.   

Figure D1 presents the results from this implementation: the ordinate represents the real roots of

equations (D6) and (D7); the abscissa is the inverse of reduced frequency.  The three black lines (solid,

dashed, dotted) represent the loci of the roots of equation (D6); the two red lines (solid, dashed)

represent the loci of the roots of equation (D7); the open blue circle indicates the intersection of the

locus of ܺோమwith the locus of ܺூభ and is the flutter solution (a single flutter mode) for the quantities

chosen.

By inverting the horizontal coordinate of the intersection point, the flutter reduced frequency is

obtained.  By using it and the vertical coordinate of the intersection point in equation (D8), the flutter

velocity is obtained.  The results are contained in Table DI.

(D8)
72

Table DI

3DOF Flutter Prediction for

Standard Case

vf kf

fps        

 

179.49    0.4476   

The problem‐specific quantities for the illustrative examples in Appendices C and D of the present paper

are the same, including the dimensional values of  and h.  As stated in Appendix C of the present

paper, alternate solution method no. 1 for the two‐degree‐of‐freedom (2DOF) flutter equation is directly

analogous to the three‐degree‐of‐freedom (3DOF) solution method.  Thus, in comparing these two

illustrative examples, one might expect their results to be similar in some ways, which is, in fact, the

case.

Comparing first the values of flutter velocity and flutter reduced frequency from the two illustrative

examples, from Tables CI and DI, it is seen that the respective values of vf and kf differ by only 3.5% and

2.7%.  These differences are attributed, obviously, to the presence of the third degree of freedom in the

3DOF solution.

Comparing next the loci of the roots in figures C1 and D1, an obvious difference is the number of loci in

each figure – three in the former and five in the latter.  However, in figure D1, if one ignores the loci of

the third “real” root,	ܺோయ, and the second “imaginary” root, ܺூమ, (those attributed to the presence of the

third degree of freedom), the remaining loci agree very well in their general character with the loci in

figure C1.  In each figure:  (1) with increasing values of the inverse of reduced frequency, the loci of the

first and second roots from the higher‐order (“real”) equation (shown in solid and dashed black)

approach each other; (2) the locus of the root from the lower‐order (“imaginary”) equation (shown in

solid red) lies between the other two loci; and (3) the locus of the “imaginary” root intersects the locus

of the second “real” root.

Results from 3DOF Solution Method Asymptotically Approach Results from 2DOF Solution Method

Three examples are presented for the 3DOF solution method.  In each example, two of the three modal

frequencies (, , h) are fixed and the third is varied from its initial value to a very high value.  In the

limit, each example corresponds to one of the three subcases discussed in the main body of the present

paper.  For example, for Subcase 1, which has as degrees of freedom  and h,  and h are fixed while
73

 is varied.  With the value of the third modal frequency at least an order of magnitude higher than the

highest value of the other two, the expectation is that the flutter velocities (vf) and flutter reduced

frequencies (kf) obtained using the 3DOF solution method will asymptotically approach those quantities

obtained using the appropriate 2DOF solution method.

Results from 2DOF solution method. – The 2DOF solution methods are solved as described in the main

body of this paper with quantities r and rr defined according to which of the three 2DOF solution

methods is employed.

Each of these examples employs the problem‐specific quantities from the NACA 496 “standard case”:

 = 0.1;  c = 0.5;  a = ‐0.4;  x = 0.2;  r

2 = 0.25;  x = 1⁄80;  r2 = 1⁄160 

The following table contains the natural frequencies and the corresponding values of frequency ratios

for the three subcases.

Table DII

Summary of Natural Frequencies and Frequency Ratios for Examples

  h   h

   rps rps rps

Subcase 1 100 x 50 1 x 1

Subcase 2 x 44.721 50 x 0.005 1

Subcase 3 100 75 x 71.111 1 x

To obtain the 2DOF flutter solutions for Subcase 1, figures 9 and 10 are used with h = 1; for Subcase 2,

figures 7 and 8 are used with  = 0.005; and for Subcase 3, figures 5 and 6 are used with  = 71.111.   

Figures D2(a) and D2(b) contain close‐up views of portions of figures 9 and 10; figures D3(a) and D3(b)

contain close‐up views of portions of figures 7 and 8; and figures D4(a) and D4(b) contain expanded

views of figures 5 and 6.  Following the convention established in the main body of the present paper,

these curves are red, indicating they were created using present calculations.  In each figure, the solid

black line indicates the target value of its particular .  Intersections of the solid black line with the red

curves produce the dashed black lines, which indicate the corresponding values of 1/kf and F, from

which kf and vf, respectively, are obtained.  (Recall that in figures 6, 8, and 10, although labelled “F” the

quantity plotted is actually the square of the defined quantity, which has ramifications on the extraction

of vf from F.)

In figures D2(a) and D2(b), there is only a single intersection of the black line with the red curves,

indicating only one flutter mode is predicted for Subcase 1; but in figures D3(a) and D3(b) and in figures
74

D4(a) and D4(b), because there are two intersections of the black line with the red curves, two flutter

modes each are predicted for Subcases 2 and 3.

The following table presents a summary of the 2DOF flutter predictions:

Table DIII

2DOF Flutter Predictions for Standard Case

Mode 1 Mode 2

vf kf vf kf

fps fps

Subcase 1 173.26 0.4355 x x

Subcase 2 19.521 2.587 120.65 0.4727

Subcase 3 14.668 8.045 234.05 0.4458

Results from 3DOF solution method. – Figures D5, D6, and D7 contain the results of many 3DOF flutter

solutions.  These figures are comprised of semi‐log plots of vf and kf as functions of the appropriate to

the subcase.  The solid blue circles represent the values of vf and kf predicted using the 3DOF solution

method, where each corresponding pair came directly from an intersection such as the one indicated by

the open blue circle in figure D1.  The dashed red lines in the figures represent the values of vf and kf

predicted using the 2DOF solution method and contained in Table DIII.   

The nominal starting value for variable frequencies , , and h in figures D5 through D7 was 100 rps.  

However, for Subcases 1 and 3, the fixed value of  is also 100 rps.  For these two subcases, the flutter

modes that ultimately became the asymptotic solutions did not develop until their respective variable

frequencies were above 100 rps:  125 rps for Subcase 1; 175 rps for Subcase 3.

As can be seen in figure D5, the 3DOF solution method predicts a single flutter mode, as did the 2DOF

solution method.  Also, with increasing values of , the circles approach the dashed lines, indicating

that, in the limit, the values of vf and kf predicted by the 3DOF solution method approach those

predicted by the 2DOF solution method, thereby confirming the initial expectation.

As can be seen in figure D6, the 3DOF solution method predicts two flutter modes, as did the 2DOF

solution method.  In each plot, with increasing values of , the circles approach the dashed lines, again

indicating that, in the limit, the values of vf and kf predicted by the 3DOF solution method approach

those predicted by the 2DOF solution method, and again confirming the initial expectation.  The
75

comparatively higher density of circles for this subcase between  = 100 rps and  = 250 rps was

necessary in order to capture the region of high curvature and the inflection in figure D6(c).

As can be seen in figure D7, the 3DOF solution method again predicts two flutter modes, as did the

2DOF solution method.  In each plot, with increasing values of h, the circles approach the dashed lines,

once again indicating that, in the limit, the values of vf and kf predicted by the 3DOF solution method

approach those predicted by the 2DOF solution method, and once again confirming the initial

expectation.

For each subcase, in the limit, the 3DOF solution method predicts the same number of flutter modes

and the same values of flutter velocity and flutter reduced frequency as did the 2DOF solution method.   

   
76

REFERENCES

1. Zeiler, Thomas A.:   Results of Theodorsen and Garrick Revisited.  Journal of Aircraft. Vol. 37, No.

5, Sept‐Oct 2000, pp 918‐920.

2. Theodorsen, Theodore:  General Theory of Aerodynamic Instability and the Mechanism of

Flutter.  NACA Report No. 496, 1940.

3. Theodorsen, Theodore; and Garrick, I. E.:  Mechanism of Flutter, a Theoretical and Experimental

Investigation of the Flutter Problem.  NACA Report No. 685, 1940.

4. Theodorsen, Theodore; and Garrick, I. E.:  Flutter Calculations in Three Degrees of Freedom.  

NACA Report No. 741, 1942.

5. Bisplinghoff, R. L.; Ashley, H.; and Halfman, R. L.:  Aeroelasticity, Addison‐Wesley‐Longman,

Reading, MA, 1955.

6. Bisplinghoff, R. L.; and Ashley, H.:  Principles of Aeroelasticity, Dover, New York, 1975.

   
77

TABLE I. ‐ TRIPPING POINTS IN NACA 496

No.     Tripping Point Identified herein

on page ‐ 

Occurs in

NACA 496 ‐ 

1 Three versions of NACA

496 dated 1934, 1940,

and 1949.

4  ‐‐ 

2 No errata for some of

the errors in 1934

version that were

corrected in 1940 and

1949 versions.

4 ‐‐ 

3 Assignment of the

same equation

numbers, (1) ‐ (9) and

(11) ‐ (19), to different

expressions.

12 pp. 12, 14, 15

4 The use of the same

symbol, Iij, to represent

a quantity and that

quantity divided by

reduced frequency.

13 pp. 12 & 14

5 The use of the same

symbol, F,  to represent

a quantity and the

square of that quantity.

16 figs 6, 8, 10 –

the square of

the quantity

‐‐  

figs 11‐14 –

the quantity

6 Quantity D is never

defined.  

18 D is plotted

in fig 12

7 Omission of values for

h and  .

23 these values

required to

produce figs

16 and 17

   
78

TABLE II. – IDENTIFIED ERRORS IN NACA 496

No. Identified Error

  

Identified herein

on page(s) ‐ 

Occurs in

NACA 496

on page ‐ 

1 The ܴ௔୦

ᇱᇱ row of Table III corresponding to

a = ‐0.4 and c = 0.5 is in error.

15 17

2 Expression ݒ ஽ൌ 	ݒோ

ଵ

భ

మ

ା௔  is incorrect.

Should be ݒ ஽ൌ 	ݒோ

ଵ

√ଵାଶ௔ .

19 16

3 As presented in figure 13, the values

comprising curve B are the square roots

of the proper values.

20 20

4 Transcription errors:

As presented in figure 14, correct values

of F were paired with non‐corresponding

values of xa – that led to incorrect curves

A, B, and C.

22 20

   
79

TABLE III. – COMPUTATIONAL SHORTCUTS IN NACA 496

No. Computational Shortcut Benefit Identified herein

on page(s) ‐ 

Occurs in

NACA 496

on page(s) ‐ 

1 Separating an equation with

complex coefficients into its

real and imaginary parts,

yielding two equations with

real coefficients.

Avoids complex arithmetic. 9, 10, 66, 67,

and 70

12 and 13

2 Eliminating the factor 1/b

from equation (10/C),

resulting in equations (12/3),

(12/6), (12/9), (12/13),

(12/16), and (12/19).

Relieves engineers of extra,

but unnecessary,

multiplications in the hand‐

calculation of these

expressions

12 12

3 Eliminating the factor 1/k

from equations (12/11)

through (12/19), resulting in

equations (14/11) through

(14/19).

Relieves engineers of extra,

but unnecessary,

multiplications in the hand‐

calculation of these

expressions

12 14
Figure 1. – Parameters of the airfoil‐aileron combination

80
Figure 2. – Constants resulting from the integration of velocity potentials.

81
Figure 3. – Elements of the Theodorsen circulation function.

82
Figure 4. – Reduced‐frequency‐dependent portions of the elements of the equations of motion.

(a) Elements from pitching moment equation.

83
Figure 4. – Reduced‐frequency‐dependent portions of the elements of the equations of motion.

(b) Elements from hinge moment equation.

84
Figure 4. – Reduced‐frequency‐dependent portions of the elements of the equations of motion.

(c) Elements from vertical force equation.

85
Figure 4. – Reduced‐frequency‐dependent portions of the elements of the equations of motion.

(d) Element           ܴ re‐drawn. ௔௛

ᇱᇱ

86
FIGURE 5. – Case 3, Torsion-aileron (, ): Standard case. Showing  against 1/k.

FIGURE 6. – Case 3, Torsion-aileron (, ): Standard case. Showing flutter factor

F against 

87
FIGURE 7. – Case 2, Aileron-deflection (h): (a) Standard case. (b), (c), (d) indicate

dependency on x. Case (d), x = -0.004, reduces to a point.

FIGURE 8. – Case 2, Aileron-deflection (h): Final curves giving flutter factor F

against  corresponding to cases shown in figure 7.

88
FIGURE 9. – Case 1, Flexure-torsion (h, ): Standard case. Showing h against 1/k.

FIGURE 10. – Case 1, Flexure-torsion (h, ): Standard case. Showing flutter factor

F against h.

89
FIGURE 11. – Case 1, Flexure-torsion (h, ): Showing dependency of F on h/. The

upper curve is experimental. Airfoil with r = 1⁄2; a = -0.4; x = 0.2; 4 = 0.01; 1/ variable.

FIGURE 12. – Case 1, Flexure-torsion (h, ): Showing dependency of F on location

of axis of rotation a. Airfoil with r = 1⁄2; x = 0.2;  = 1⁄4; 1/ = 1⁄6; a variable.

90
Figure 13A. – Curve A

FIGURE 13. – Case 1, Flexure-torsion (h, ): Showing dependency of F on the radius

of gyration r = r. 

A, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1⁄6; r variable.

B, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1.00; r variable.

91
Figure 13B. – Curve B

(a) Comparison of present calculation and Curve B of NACA 496

FIGURE 13. – Case 1, Flexure-torsion (h, ): Showing dependency of F on the radius

of gyration r = r. 

A, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1⁄6; r variable.

B, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1.00; r variable.

92
Figure 13B. – Curve B

(b) Comparison of present calculation and square of Curve B of NACA 496

FIGURE 13. – Case 1, Flexure-torsion (h, ): Showing dependency of F on the radius

of gyration r = r. 

A, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1⁄6; r variable.

B, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1.00; r variable.

93
Figure 13B. – Curve B

(c) Comparison of present calculation, independent calculation, and square of Curve B of NACA 496

FIGURE 13. – Case 1, Flexure-torsion (h, ): Showing dependency of F on the radius

of gyration r = r. 

A, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1⁄6; r variable.

B, airfoil with a = -0.4;  = 1⁄4; x = 0.2; 1/ = 1.00; r variable.

94
Figure 14A. – Curve A

(a) Comparison of present calculation and Curve A of NACA 496

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

95
(b) Comparison of present calculation and rotated Curve A of NACA 496

Figure 14A. – Curve A

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

96
(c) Comparison of present calculation, independent calculation, and rotated Curve A of NACA 496

Figure 14A. – Curve A

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

97
Figure 14B. – Curve B

(a) Comparison of present calculation and Curve B of NACA 496

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

98
Figure 14B. – Curve B

(b) Comparison of present calculation and rotated Curve B of NACA 496

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

99
(c) Comparison of present calculation, independent calculation, and rotated Curve B of NACA 496

Figure 14B. – Curve B

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

100
Figure 14C. – Curve C

(a) Comparison of present calculation and Curve C of NACA 496

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

101
Figure 14C. – Curve C

(b) Comparison of present calculation and rotated Curve C of NACA 496

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

102
(c) Comparison of present calculation, independent calculation, and rotated Curve C of NACA 496

Figure 14C. – Curve C

FIGURE 14. – Case 1, Flexure-torsion (h, ): Showing dependency of F on x, the

location of the center of gravity.

A, airfoil with r = 1⁄2; a = -0.4;  = 1/400; 1/ = 1⁄6; x variable.

B, airfoil with r = 1⁄2; a = -0.4;  = 1⁄4; 1/ = 1⁄6; x variable.

C, airfoil with r = 1⁄2; a = -0.4;  = 1/100; 1/ = 1; x variable.

103
FIGURE 15. – Case 1. Wing A. Theoretical and experimental curves giving flutter

velocity v against frequency ratio h/. Deflection-torsion.

104
-- with assumed value of h = 5.8x2 --

-- with assumed value of  = 4.4x2 --

FIGURE 16. – Case 2. Wing B. Theoretical and experimental curves giving flutter

velocity v against frequency ratio /h. Aileron-deflection (, h).

FIGURE 17. – Case 3. Theoretical curve giving flutter velocity against the frequency ratio /. The experimental unstable are is indefinite due to the importance of internal friction at very small velocities. Torsion-aileron (, ).

105
Figure C1 – Illustrative example for Alternate Solution Method No. 1.

106
Figure C3 – Illustrative example for Alternate Solution Method No. 2

– Reduced frequencies as a function of velocity.

Figure C2 – Illustrative example for Alternate Solution Method No. 2

– Z1 and Z2 as functions of reduced frequency at a velocity of 160 fps.

107
Figure C4 – Illustrative example for Alternate Solution Method No. 3.

108
Figure D1 – Illustrative example for Three‐Degree‐of‐Freedom Solution

Method; “standard case” quantities chosen, with  = 100,  = 125, h = 50.

109
Figure D4 – Results from 2DOF Solution Method, Subcase 3;

 = 100,  = 75, r

2 = 0.25, r2 = 1/160, x = 1/80,  = 71.111.

(a) Recreation of Figure 5 (b) Recreation of Figure 6

Figure D3 – Results from 2DOF Solution Method, Subcase 2;

 = 44.721, h = 50, r

2 = 1/160, x = 1/80,  = 0.005.

(a) Recreation of Figure 7 (b) Recreation of Figure 8

Figure D2 – Results from 2DOF Solution Method, Subcase 1;

 = 100, h = 50, r

2 = 0.25, h = 1.

(a) Recreation of Figure 9 (b) Recreation of Figure 10

110
Figure D5 – Results from 3DOF Solution Method, “Subcase 1”;

 = 100, h = 50,  variable, “standard case”.

(a) Flutter velocity vs  (b) Flutter reduced frequency vs 

Figure D6 – Results from 3DOF Solution Method, “Subcase 2”;

 = 44.721, h = 50, r

2 = 1/160, x = 1/80,  variable.

(a) Mode 1 flutter velocity vs  (b) Mode 1 flutter reduced frequency vs 

(c) Mode 2 flutter velocity vs  (d) Mode 2 flutter reduced frequency vs 

111
Figure D7 – Results from 3DOF Solution Method, “Subcase 3”;

 = 100,  = 75, r

2 = 0.25, r2 = 1/160, x = 1/80, h variable.

(a) Mode 1 flutter velocity vs h (b) Mode 1 flutter reduced frequency vs h

(c) Mode 2 flutter velocity vs h (d) Mode 2 flutter reduced frequency vs h

112
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Publication

 4. TITLE AND SUBTITLE

Re-Computation of Numerical Results Contained in NACA Report No. 496

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Perry, Boyd, III

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20554

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 08

Availability: STI Program (757) 864-9658

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

An extensive examination of NACA Report No. 496 (NACA 496), “General Theory of Aerodynamic Instability and the Mechanism of Flutter,” by Theodore 

Theodorsen, is described. The examination included checking equations and solution methods and re-computing interim quantities and all numerical examples 

in NACA 496. The checks revealed that NACA 496 contains computational shortcuts (time- and effort-saving devices for engineers of the time) and clever 

artifices (employed in its solution methods), but, unfortunately, also contains numerous tripping points (aspects of NACA 496 that have the potential to cause 

confusion) and some errors. The re-computations were performed employing the methods and procedures described in NACA 496, but using modern 

computational tools. With some exceptions, the magnitudes and trends of the original results were in fair-to-very-good agreement with the re-computed results. 

The exceptions included what are speculated to be computational errors in the original in some instances and transcription errors in the original in others. 

Independent flutter calculations were performed and, in all cases, including those where the original and re-computed results differed significantly, were in 

excellent agreement with the re-computed results. Appendix A contains NACA 496; Appendix B contains a Matlab ® program that performs the 

re-computation of results; Appendix C presents three alternate solution methods, with examples, for the two-degree-offreedom solution method of NACA 496; 

Appendix D contains the three-degree-of-freedom solution method (outlined in NACA 496 but never implemented), with examples. 

15. SUBJECT TERMS

Aeroelastic stability; Flutter; NACA Report No. 496; Theodore Theodorsen

18. NUMBER

 OF 

 PAGES

117

19b. TELEPHONE NUMBER (Include area code)

(757) 864-9658

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

432938.11.01.07.43.40.08 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TP-2015-218765

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 06 - 2015