/I* , -- 

A "ERICAL MODEL OF UNSTEADY, 

SUBSONIC AEROELASTIC REHAVIOR 

BY 

Thomas W. Strganac 

Virginia Polytechnic Institute and State University 

in Partial Fulfillment of the Requirements for the Degree of 

DOCTOR OF PHILOSOPHY 

in 

Engineering Mechanics 

APPROVED : 

Dean T. Mook, Chairman 

Daniel Frederick Leonard Meirovitch 

Scott L. Hendricks William F. Hunter 

August, 1987 

Blacksburg, Virginia 

(NASA-TH-1Q1126) A ?irJfiEI?ICAL NODEL OF N a 8- 2 28 62 

UNSTEADY, SUBS3NIC AEROELASTIC REHRVIOR Ph, D 

Thesis (Virqinia Polytechnic Inst. akd 

State aniv.) 123 p CSCL ula Unclas 

(;3/02 0 142549 
A NUMERICAL MODEL OF UNSTEADY, 

SUBSONIC AEROELASTIC BEHAVIOR 

BY 

Thomas W. Strganac 

ABSTRACT 

A method for predicting unsteady, subsonic aeroelastic responses 

has been developed. The technique accounts for aerodynamic nonlinearities associated with angles of attack, vortex-dominated flow, static 

deformations, and unsteady behavior. The angle of attack is limited 

only by the occurrence of stall or vortex bursting near the wing. The 

fluid and the wing together are treated as a single dynamical system, 

and the equations of motion for the structure and flowfield are integrated simultaneously and interactively in the time domain. The 

method employs an iterative scheme based on a predictor-corrector 

technique. The aerodynamic loads are computed by the general unsteady 

vortex-lattice method and are determined simultaneously with the motion of the wing. Because the unsteady vortex-lattice method predicts 

the wake as part of the solution, the history of the motion is taken 

into account; hysteresis is predicted. Two models are used to demonstrate the technique: a rigid wing on an elastic support experiencing 

plunge and pitch about the elastic axis, and an elastic wing rigidly 

supported at the root chord experiencing spanwise bending and twisting. The method can be readily extended to account for structural 

nonlinearities and/or substitute aerodynamic load models. The time 

domain solution coupled with the unsteady vortex-lattice method 

provides the capability of graphically depicting wing and wake motion. 
ACKNOWLEDGMENTS 

My experience at Virginia Tech was genuinely worthwhile because 

of the individuals - students, staff, and faculty - who exemplify the 

university. I wish to thank the members of my graduate committee - 

Professors Frederick, Meirovitch, and Hendricks - who assisted me in 

my graduate program and provided an excellent educational opportunity 

in the classroom. I am grateful to Mrs. Sally Shrader who typed 

several conference papers on my behalf and became my liaison with the 

university in my absence. 

I am especially indebted to my mentor, Professor Mook, who assured my education went well beyond the classroom and was not limited 

by the research contained within this dissertation. His patience and 

attitude was certainly tested through our collaboration. Yet I always 

came away from each of our discussions with an improved understanding 

of our research and an appreciation of Professor Mook's sincere 

dedication to his profession. 

I am grateful that NASA allowed me the opportunity to participate 

in the Graduate Study Program. Yet, this study could not have been 

realized without the continuous support of the staff of the Langley 

Research Center. Individuals throughout the organization such as the 

members of the Employee Development Branch, Visual and Printing 

Services Branch, Computer Applications Branch, and Flight Software and 

Graphics Branch have all contributed to a successful program. I am 

grateful to Dr. William F. Hunter who supported my interest in this 

research and served on my graduate committee. I am indebted to 

Mrs. Maria Vallas Mitchum who provided assistance in programming and 

iii 
graphics. And I wish to thank the members of the Technical Editing 

Branch and, in particular, Miss Patricia Gottschall €or the typing of 

this manuscript. 

Finally, I am thankful to my wife, Kathy, who supports my 

interests at her sacrifice and to our children, Christopher and Kasey, 

who provide the purpose. 

iv 
TABLE OF CONTENTS 

ACKNOWLEDGEMENTS ................................................ iii 

TABLE OF CONTENTS ............................................... V 

NOMENCLATURE .................................................... X 

CHAPTER I . INTRODUCTION ........................................ 1 

1.1 General 1 

1.2 The Aeroelastic Phenomenon ............................... 2 

1.3 Literature Review ........................................ 7 

1.3.1 Aeroelastic models ................................ 7 

.................................................. 

1.3.2 Aerodynamic models ................................ 11 

1.3.3 Structural models ................................. 14 

1.4 The Present Method ....................................... 15 

CHAPTER I1 . EQUATIONS OF MOTION FOR THE RIGID WING ............. 17 

2.1 Description of the System ................................ 17 

2.2 Development of the Equations ............................. 19 

2.3 Nondimensionalized Form of the Governing Equations ....... 23 

CHAPTER I11 . EQUATIONS OF MOTION FOR THE ELASTIC WING .......... 25 

3.1 Description of the System ................................ 25 

3.2 Development of the Equations ............................. 25 

3.3 Static and Dynamic Equations ............................. 29 

3.3.1 Equations governing the static contribution ....... 30 

3.3.2 Equations governing the dynamic contribution ...... 33 

CHAPTER IV . AERODYNAMIC LOADS. THE UNSTEADY 

VORTEX-LATTICE METHOD .............................. 4.1 Overview ................................................. 4.2 A Description of the Method .............................. 4.2.1 The wing representation ........................... 4.2.2 Kinematic considerations .......................... 4.2.3 The Biot-Savart Law ............................... 4.2.4 The boundary conditions ........................... 4.2.5 The calculation of the circulation ................ 4.2.6 The unsteady Kutta condition ...................... 4.3 The Impulsive Start of the WLM .......................... 4.4 Graphical Representations of the Flowfield ............... 4.5 The Aerodynamic Loading .................................. 

36 

36 

37 

37 

39 

40 

41 

42 

43 

43 

44 

44 

V 
CHAPTER V. THE NUMERICAL SOLUTION SCHEME ....................... 5.1 Overview ................................................. 5.2 The First-Order Differential Equations ................... 5.3 The Basic Predictor-Corrector Method ..................... 5.4 Starting the Predictor-Corrector Method .................. 5.5 Integration Method for the Rigid-Wing System ............. 5.5.1 First-order convection theory ..................... 5.5.2 Second-order convection theory .................... 5.6.1 Static solution ................................... 5.6.2 Dynamic solution .................................. 5.6 Solution of the Equations for the Elastic Wing ........... 

CHAPTER VI . APPLICATION OF THE MODEL TO 

SPECIFIC EXAMPLES .................................. 6.1 Introduction ............................................. 6.2 An Example of Aeroelastic Behavior €or 

the Rigid Wing ........................................... 6.2.1 Subcritical response without 

structural damping ................................ 6.2.2 Supercritical response to a small 

initial disturbance ................................ 6.2.3 Supercritical response to a large 

initial disturbance ............................... 6.2.4 Subcritical response with 

structural damping ................................ 6.2.5 Comparison with theoretical results ............... An Example of Aeroelastic Behavior for 

the Elastic wing..............................^^^^^^^^^^^^. 

6.3.1 Subcritical response of the elastic wing .......... 6.3.2 Critical response of the elastic wing ............. 6.3.4 Comparison with empirical results ................. 6.3.5 Unsteady aerodynamic loads ........................ 6.4 The Identification of Flutter ............................ 6.4.1 Parametric identification procedures .............. 6.4.2 Identification using the system energy ............ 

6.3 

6.3.3 Supercritical response of the elastic wing ........ 

6.5 Wing and Wake Graphics ................................... CHAPTER VI10 CONCLUDING REMARKS ................................ REFERENCES ...................................................... VITA ............................................................ 

49 

49 

49 

51 

53 

56 

56 

58 

58 

58 

60 

64 

64 

64 

66 

68 

49 

72 

74 

76 

78 

81 

53 

85 

86 

91 

92 

94 

97 

101 

102 

110 

vi 
LIST OF FIGURES 

Figure 1.1. The Three-Ring Diagram. 4 

Figure 1.2. Flutter Response for a Two Degree of Freedom Example. 6 

Figure 1.3. Determination of Flutter for Frequency Domain 

Formulations, V-g and Frequency Diagrams. 9 

Figure 2.1. The Rigid-Wing Model. 18 

Figure 3.1. The Elastic-Wing Model. 26 

Figure 4.1. The Unsteady Vortex-Lattice Method. 38 

Figure 4.2. Wireframe of Wing and Wake Lattice. 45 

Figure 4.3. Vortex-Lattice Representations. 46 

Figure 4.4. Wing and Wake Vorticity for Wings of Different 

Aspect Ratios. 47 

Figure 5.1. Flowchart of the Integration Scheme Using LowerOrder Convection Theory, Rigid-Wing Formulation. 57 

Figure 5.2. Flowchart of the Integration Scheme Using SecondOrder Convection Theory, Rigid-Wing Formulation. 59 

Figure 5.3. Flowchart for the Determination of the Elastic 

Wing Deformations Due to Static Loads. 61 

Figure 5.4. The Integration Process for the Equations Governing 

the Dynamic Solution of the Elastic Wing. 62 

Figure 6.1. Response of the Rigid Wing at a Subcritical 

Dynamic Pressure. 67 

Figure 6.2. Response of the Rigid Wing at a Critical Dynamic 

Pressure - Small Initial Conditions. 70 

Figure 6.3. Response of the Rigid Wing at a Critical Dynamic 

Pressure - Large Initial Conditions. 71 

Figure 6.4. Phase Plane Diagrams at the Critical Dynamic 

Pressures - Small and Large Initial Conditions. 73 

Figure 6.5. Response of the Rigid Wing at a Critical Dynamic 

Pressure With Structural Damping Present. 75 

Figure 6.6. Response of the Elastic Wing at a Subcritical 

Dynamic Pressure. 80 

vii 
Figure 6.7. Response of the Elastic Wing at a Critical 

Dynamic Pressure. 

Figure 6.8. Response of the Elastic Wing at a Supercritical 

Dynamic Pressure. 

Figure 6.9. The Pressure Distribution of the Elastic Wing 

at Three Different Times. 

Figure 6.10. The Predicted Pressure Distribution for 

Four Selected Elements. 

Figure 6.11. Sectional Aerodynamic Coefficients Predicted 

by the Model. 

Figure 6.12. The Computed Energy of the Rigid Wing for a 

Subcritical and Critical Case. 

Figure 6.13. The Strength of Vorticity for the Wing at 

Several Time Steps. 

Figure 6.14. The Predicted Shape and Strength of Vorticity 

for the Wing and Wake at Several Time Steps. 

82 

84 

87 

89 

90 

96 

99 

100 

viii 
LIST OF TABLES 

Table 6.1. Properties of the Rigid Wing Example 65 

Table 6.2. Properties of the Elastic Wing Example 77 

Table 6.3. Natural Modes Used to Represent the Elastic Wing 79 

ix 
NOMENCLATURE 

A 

AR 

(AI 

C 

‘a 

‘m 

Cn 

C 

‘a 

Y 

E1 

C 

FX 

FY 

GJ 

g 

Icg 

‘e 

[I1 

K 

Ka 

5 

[K I 

ka 

wing reference area 

aspect ratio 

generalized aerodynamic loading vector 

nondimensional aerodynamic constant = p tCAI2 m 

nondimensional rotational damping = ca/m UmIlc 

aerodynamic pitch-moment coefficient 

aerodynamic normal-force coefficient 

nondimensional translational damping = c 

chord 

support rotational damping 

support translational damping 

distributed bending stiffness 

horizontal component of the reaction of the support structure 

vertical component of the reaction of the support structure 

distributed torsional stiffness 

acceleration due to gravity 

mass moment of inertia referenced to the center of gravity 

mass moment of inertia referenced to the elastic axis 

identity matrix 

11 /m UYC 

nondimensional 

nondimensional 

nondimensional 

2 inertia = Ie/m kc 

2 

rotational stiffness = ka/m U m 

2 2 translational stiffness = k II /m UoD YC 

stiffness matrix 

support rotational stiffness 

X 
kY 

IlC 

M 

[MI 

Mr 

MS 

m 

N 

n 

- 

X 

Y 

a 

support translational stiffness 

characteristic length = c/n 

distributed aerodynamic pitch moment 

mass matrix 

aerodynamic pitch moment at the reference position 

moment due to support structure 

mass or distributed mass 

distributed aerodynamic normal force 

number of chordwise (bound) lattice elements or number of 

first-order equations 

generalized coordinates 

distance between center of gravity (cg) and elastic axis 

(ea), positive for cg aft of ea 

static loading vector 

freestream velocity 

displacement of the elastic axis 

local distance between the sectional elastic axis and center 

of gravity 

distance between center of gravity (cg) and the aerodynamic 

reference, positive for cg aft of reference 

C 

nondimensional plunge degree of freedom = y/ll 

nondimensional pitch degree of freedom, or rotation about the 

elastic axis 

time derivative of Y or the generalized velocities (= (;I 

time derivative of a 

density of fluid 

bending or torsional free-vibration mode 

static angle-of-attack 

xi 
( )' first derivative with respect to length 

( )*I second derivative with respect to length 

(*I first derivative with respect to time (or nondimensional 

time 

(**> second derivative with respect to time (or nondimensional 

time 1 

xii 
CHAPTER I 

INTRODUCTION 

1.1 General 

Aeroelastic instabilities, which may be catastrophic, occur as a 

result of coupling between aerodynamic forces, structural forces, and 

mass. Therefore, the aeroelastic design of aircraft must address the 

structure and the aerodynamic forces within the flight environment. 

This concern includes the design of wings or aircraft for both flight 

and wind-tunnel testing. 

A method to predict the aeroelastic behavior of a wing is described. The objective is to develop an improvement in aeroelastic 

analyses by accounting for aerodynamic nonlinearities associated with 

angles of attack, static deformations, vorticity-dominated flow, and 

unsteady behavior. This method generates a realistic simulation of 

aeroelastic response or instability by predicting both the motion of 

the wing and the motion of the fluid simultaneously. In other words, 

the wing and fluid are treated as a single dynamical system, and the 

equations of motion for the structure and flowfield are integrated 

simultaneously and interactively in the time domain. 

The equations of motion are developed for classical two-degreeof-freedom wing motion; that is, a rigid wing which is allowed to 

plunge and pitch about the elastic axis. In addition, the equations 

of motion are developed for an elastic wing with bending and twisting. 

The unsteady vortex lattice method (UVLM) is used to predict the 

aerodynamic loads. The technique accounts for the aerodynamic 

nonlinearities. 

1 
2 

A complication does exist with a time domain approach: the aerodynamic loads cannot be predicted unless the motion of the wing is 

known, and the motion of the wing cannot be predicted unless the aerodynamic loads are known. An iterative scheme that accounts for the 

interaction between the aerodynamic loads and wing motion was developed based upon a predictor-corrector approach. An advantage of the 

method is that it provides an excellent opportunity to animate the 

motions in the flowfield and wing slnce the governing equations are 

solved in the time domain. 

The technique described herein provides a somewhat different approach to model flutter for a wing of arbitrary planform. 

tions are developed about arbitrary static angles of attack, providing 

a new capability to study the associated nonlinear effects. Through 

the use of the UVLM one is able to treat low-aspect-ratio wings and 

account for the wake rollup at the wing tips. The UVLM models the 

wake where the history of the motion resides. 

The equa1.2 The Aeroelastic Phenomenon 

Many different types of aeroelastic instabilities exist. Static 

instabilities include aeroelastic divergence in which the elastic restoring forces of the wing are exceeded by the aerodynamic forces. 

Dynamic instabilities such as flutter are caused by an exchange of 

energy between the aircraft and surrounding air. General aeroelastic 

phenomena are described by Fung (1955), Bisplinghoff, Ashley and 

Halfman (19551, and Dowel1 (1980) among others. A historical perspective of aeroelastic research is given by Collar (1978). 
3 

A complete 

inertial forces 

model of the aerodynamic forces, elastic forces, and 

is requited to describe the unsteady aeroelastic nature of an aircraft. As Figure 1.1 depicts, this interaction among 

these three types of force provides the conditions necessary for 

potential static and dynamic instabilities. 

elastic problem can best be illustrated by this three-ring diagram as 

suggested by Yates (1971). Variations to this diagram are presented 

by Bisplinghoff et al. (19551, Collar (1978), and Dowel1 (1980). 

The nature of the aeroAerodynamic, elastic, and inertial forces are shown. The common 

ground between the inertial forces and elastic forces represents the 

free-vibration problem, the common ground between the aerodynamic 

forces and elastic forces represents static aeroelastic phenomena 

which include instabilities such as divergence, and the common ground 

between the aerodynamic forces and inertial forces represents dynamic 

stability investigations of "rigid" aircraft. The interaction of all 

three forces is present in the aeroelastic phenomenon known as 

flutter. 

Flutter has been defined in many ways. Among the first to examine the nature of flutter, Theodorsen (1940) described flutter and 

developed a general theory which is still used today. 

defines flutter as "a self-excited oscillatory aerodynamic instability" and states that flutter may be considered "an oscillatory instability where one degree of freedom is driven at resonance by a 

second degree of freedom, both at the same frequency". Hancock, 

Wright, and Simpson (1985) describe flutter as "a complex phenomenon 

where, in the classical sense of the term, two or more structural 

Pines (1958) 
4 

Figure 1.1. The Three-Ring Diagram. 
5 

normal modes are coupled and excited through time dependent aerodynamic loads". Dowel1 (1980) also describes the flutter phenomenon 

and suggests that many types of flutter exist; these include: 

"'Coalescence' or 'merging frequency' flutter", '"Single-degree-offreedom' flutter", and "'Divergence' or 'zero frequency' flutter". 

An illustration of the flutter phenomenon is presented in Figure 1.2. A wing is initially moving through a fluid at rest. The 

wing is subjected to a disturbance. In the left-hand side of the 

figure the dynamic pressure is below the critical condition for flutter. The bending and torsion of the wing occur at separate, unique 

frequencies which are close to the natural frequencies for the free 

vibration modes. The motion decays following the initial disturbance 

as there is not enough energy being extracted from the freestream to 

maintain the motion. In fact, the freestream is extracting energy 

from the wing and damping the motion. There is aerodynamic damping 

here. The system is being damped by aerodynamic effects only. In the 

right-hand side of the figure the dynamic pressure is greater than the 

critical pressure. Not only is the motion seen to grow in this case 

but bending and torsion occur at a common frequency - a characteristic 

associated with flutter. In the lower portion of the figure the shift 

in system frequencies due to the dynamic pressure is illustrated. The 

frequencies are the natural ones when the dynamic pressure is zero. A 

coalescence occurs at a critical dynamic pressure, the system is now 

in flutter. 
cr qcq cr q’q 

c b 

Figure 1.2. Flutter Response for a Two Degree of Freedom Example. 
7 

1.3 Literature Review 

The purpose of this research is to provide an improvement in 

aeroelastic analysis capabilities. 

UVLM to the aeroelastic formulation new capabilities are introduced. 

As mentioned previously these include the modeling of aerodynamic 

nonlinearities and new solution strategies. These enhancements, in 

turn, open the door for further improvements primarily in structural 

nonlinearities. However, the emphasis in this research is placed on 

the aerodynamic model and the solution strategy; the literature review 

will focus on these areas. 

Through the application of the 

1.3.1 Aeroelastic models 

Aeroelastic research is very broad indeed. Experimental and 

theoretical research has been performed in structures, fluids, and 

dynamics as well as solution schemes associated with aeroelasticity. 

Desmarais and Bennett (1978) described an approach for flutter 

analysis which is representative of current numerical strategies. 

In their method the generalized mass and mode shapes are developed 

through an external structural analysis package. 

to describe the shape of the wing which in turn is used to compute an 

aerodynamic influence matrix. A subsonic kernal function approach is 

These modes are used 

used to model the aerodynamics. An eigenvalue formulation is used 

with the governing equations. Hence, the solution is performed in the 

frequency domain, the governing equations are cast in a linear form, 

and a simple harmonic form of the solution is assumed. As we shall 

see, the method developed in this research is not limited by these 
8 

assumptions. 

upon the velocity, density, and Mach number, are used to determine the 

critical conditions necessary for flutter. 

The roots of this eigensolution, which are dependent 

The method commonly used to determine the flutter conditions of a 

wing using a frequency domain solution is illustrated in Figure 1.3. 

The method, called the V-g (for Velocity-incremental damping) method, 

tracks the roots of the modes used for the solution for increasing 

values of velocity. A fictitious damping (the root) is determined at 

each velocity which is required to satisfy the eigenvalue problem. In 

the upper part of the figure we plot the roots for a 3 mode system. 

The upper plane represents an unstable region; hence, when a root 

enters this region the system is unstable. At only this "flutter 

crossing" is the solution meaningful, yet no physical interpretation 

of the motion is available. In the lower part of the figure the 

frequencies for these modes are computed as the velocity increases. 

Hassig (1971) and Fung (1955) described some alternate solution 

schemes for the frequency domain solutions. Fung also presented many 

examples, some of which will be referred to in this dissertation. 

Goland and Luke (1949) outlined a method to describe the aeroelastic 

character of an elastic wing. 

Garrick (1969) presented several classical papers on the flutter 

phenomenon, flutter determination, aerodynamic theories, and other 

associated subjects. Loring (1941) introduced generalized coordinates 

through the use of the natural modes to solve the flutter problem. 

His work developed a systemized approach for the structural model. 
9 

UNSTABLE I 

MODE 1 

constant 

Density 

Mach No 

\ 

MODE 1 

FLUTTER SPEED INDEX 

Figure 1.3. Determination of Flutter for Frequency Domain 

Formulations, V-g and Frequency Diagrams. 
10 

Dowel1 (1980) described both frequency and time domain solutions 

of the aeroelastic equations. Frequency domain solutions are the most 

widely used since these require less resources and a more convenient 

solution scheme. However, a frequency domain solution does not provide a physical description of the motion as do time domain solutions. 

Time domain solutions, which do not assume a form of the solution 

(i.e. harmonic motion) as do the frequency domain schemes, provide the 

more general and physical approach. 

A more complete formulation of the aeroelastic equations was 

described by Bisplinghoff (1975) and Fung (1955). Descriptions of 

aeroelastic phenomena were presented by Theodorsen (19401, Biot and 

Arnold (19481, Pines (1958), and Hancock et al. (1985). These later 

two papers presented physical models of flutter. Pines, in particular, examined the governing equations in a general sense and 

established relationships for which flutter can exist. 

An aeroelastic model, which is similar in nature to the one presented in this dissertation, is described by Devers (1972). Devers 

developed a time domain solution of the governing equations, he too 

used a vortex-lattice method as an aerodynamic model. However, simplifications in his model prevented the unsteady nature of trailing 

edge and wing-tip vortex effects to be considered. In addition, his 

integration scheme assumed a form of the solution and, as a consequence, the simulation reflected this assumption. 
11 

1.3.2 Aerodynamic models 

Currently, the major effort of aeroelastic research is to describe the unsteady aerodynamic loads, and, in particular, the flow in 

the transonic flow regime as most recently described by Batina et al. 

(1987). Edwards (1986) presented an overview of aerodynamics associated with transonic aeroelasticity. 

Bisplinghoff (1975) described many different aerodynamic operators which are available for all flow regimes. Yates (1985) described 

the various approaches available for aerodynamic models, which include 

doublet-lattice, kernal function, velocity potential and acceleration 

(pressure) potential formulations. 

The research will concentrate in the subsonic regime. Theodorsen 

(1940) introduced a two-dimensional functional form of the subsonic 

aerodynamic operator which is still popular today. Desmarais and 

Bennett (1978) implemented a subsonic kernal function indicative of 

more current strategies. However, most methods are quasi-steady, linear (i.e. small angle) theories. A new approach to the subsonic aerodynamic model for aeroelastic analysis is introduced. The unsteady 

vortex-lattice method (UVLM) provides the opportunity to address aerodynamic nonlinear effects associated with unsteady flow aspect ratios, 

static deformations, and the angles of attack. Hence, nonlinear 

aerodynamic effects may be studied as they apply to the aeroelastic 

phenomenon. 

Experimental investigations have demonstrated nonlinear effects. 

Yates and Bennett (1971), Farmer (19821, and Yates, Wynne, and Farmer 

(1982) addressed the effect angle of attack has on flutter boundaries. 
12 

Hsu (1957) addressed the effect aspect ratio has on flutter 

boundaries. 

The development of the WLM will be briefly reviewed. The reader 

is referred to the work of Kandil (19741, Maddox (19731, and most recently, Konstadinopoulos (1981) for more complete descriptions. 

The unsteady vortex-lattice method, as we use it in this research, 

originates from the investigation of Belotserkovskiy (1966). This 

original work could treat arbitrary, wing planforms and deformations 

but could not model the geometry of the wake which thereby limited it 

to small angles of attack. Belotserkovskiy (1968) then developed a 

nonlinear method to calculate the aerodynamic loads on wings with 

wing-tip separation in steady flows. Later, Belotserkovskiy and Nisht 

(1974) presented a method for the treating of rather general planforms 

in nonlinear unsteady flow. As part of this research, they determined 

the shape of the wake convecting from the wing tips and the trailing 

edge . 

Maddox (1973) and Mook and Maddox (1974) used this vortex-lattice 

method and considered leading edge separation. However, this method 

did not account for force-free wing-tip and trailing-edge vortex 

sheets. Kandil (1974) and Kandil, Mook, and Nayfeh (1974) refined 

this approach and their results were in very good agreement with the 

experiments. The remarkable accuracy of vortex-lattice methods was 

discussed by James (1971). 

Atta (1976) and Atta, Kandil, Mook, and Nayfeh (1976, 1977) extended the method to treat unsteady flows past rectangular wings with 

sharp-edge separation. The problem was posed in an inertial frame 
13 

which proved to be awkward. Thrasher, Mook, Kandil, and Nayfeh (1977) 

and Thrasher (1979) posed the problem in terms of a body-fixed frame. 

The approach treated rectangular wings executing arbitrary maneuvers 

as long as separation occurred along the wing's sharp edges and vortex 

bursting did not occur at or near the wing's surface. 

Thrasher (1979) coupled the aerodynamic method with a predictorcorrector method to predict the flowfield, loads, and motion of a 

hinged rectangular wing due to a prescribed motion of a flap. This 

development is of particular interest since it had application to 

aeroelastic type problems. Kandil, Atta, and Nayfeh (1977) and Atta 

(1978) refined the approach to treat delta wings. Further, they 

developed a higher-order convection process for the wake, but this 

process required iteration. 

Nayfeh et al. (1979) modified the method to treat small, harmonic 

oscillations around an arbitrary angle of attack. The general 

unsteady method was described by Nayfeh, Mook, and Yen (19791, 

Konstadinopoulos (1981) and Konstadinopoulos, Mook, and Nayfeh (1981 1. 

Most recently the general method was described by Konstadinopoulos, 

Thrasher, Mook, Nayfeh, and Watson (1985). 

Recently, Konstadinopoulos (1984) and Konstadinopoulos, Mook, and 

Nayfeh (1985) developed a numerical method for simulating subsonic 

wing rock. More recently, Elzebda (1986) described a model for simulating two-degree-of-freedom wing rock. Elzebda, Mook, and Nayfeh 

(1985) described the ability of the general method to model the aerodynamics for close-coupled lifting surfaces. 
14 

1.3.3 Structural models 

In this research the wing structure is represented by two different models for demonstration purposes - a rigid wing which can pitch 

and plunge about the elastic axis and an elastic wing. Currently, for 

the model of the elastic wing, the elastic axis of the wing is represented by a beam which is allowed to bend and twist. No structural 

nonlinearities are taken into account. The two-degree-of-freedom 

model is addressed in the texts on aeroelasticity we have already 

cited. Several texts (for example, see Rivello (1969)) describe the 

equations of a beam which we have implemented. 

The numerical model can be readily extended to include nonlinear 

effects of beam flexure and torsion. Several references consider 

nonlinear models, these include the work of Crespo da Silva and Hodges 

(19861, Hodges and Dowell (1974), and Kaza and Kvaternik (1977). 

These references address the more difficult kinematics and large 

deformations associated with rotorcraft applications. 

More complicated models of the wing are typically handled by finite element methods. For example, the method described by Desmarais 

and Bennett (1978) implements the structural modes generated by any 

experimental or analytical method. Hence, the description of structures is dependent on the capability of these alternate sources. Yet, 

an advantage in the aeroelastic design of structures is gained through 

the inclusion of more elaborate structural models. Ashley et al. 

(1980), Bendiksen and White (19861, and Weisshaar (1987) describe the 

advantages of aeroelastic tailoring of structures. 
15 

1.4 The Present Method 

A method for predicting unsteady, subsonic aeroelastic responses 

is described. 

the two-degree-of-freedom problem (Strganac and Mook (1986)). The 

introduction of the UVLM to the prediction of flutter provides the 

opportunity to model aerodynamic nonlinearities. Later, we extended 

this effort to include the elastic wing model (Strganac and Mook 

(1987)). Most recently, we demonstrated the simulation of flutter by 

animating the wing and flowfield as predicted by our model (Strganac, 

Mook, and Mitchum (1987)). These efforts are further described here. 

Previously, we described the application of the UVLM to 

In Chapter 11, the equations of motion for a rigid wing on an 

elastic foundation are developed. Structural damping and static 

deflections are included. The nondimensionalization concept is also 

introduced. 

In Chapter 111, the equations of motion for an elastic wing on a 

fixed support are developed. The equations account for static deformations. Mass and stiffness matrices for the wing are developed. 

Coupling exists between bending and torsion. Structural nonlinearities are not addressed but the general formulation allows for this 

feature. 

In Chapter IV, the unsteady vortex-lattice method is described. 

The computation of the wake geometry is illustrated. Comparisons of 

the computed pressure distribution and the sensitivity to aspect ratio 

are shown. 

In Chapter V, the integration schemes which account for the 

structural dynamic and aerodynamic interaction are developed. The 
16 

methods are based upon a predictor-corrector technique, with appropriate modifications to account for aerodynamic/dynamic coupling. 

convection schemes which are used to accurately generate the wake in 

the unsteady model are discussed. 

formulations are individually described as each present unique 

requirements. 

The 

The rigid wing and elastic wing 

In Chapter VI, several examples of computed results are presented. Wing response, aerodynamic loads, and pressure distributions 

are shown. In addition, phase planes and aerodynamic load histories 

are shown which further describe the physical nature of aeroelastic 

phenomena. Finally, the ability to take advantage of time domain 

solutions is graphically demonstrated. The wing and wake shape and 

vorticity distribution are shown at several time steps. These 

sequences are indicative of the animation capability which exists 

with the present formulation. 
CHAPTER I1 

EQUATIONS OF MOTION FOR THE RIGID WING 

2.1 Description of the System 

In this chapter the equations of motion are developed for a rigid 

wing on an elastic support limited to plunge and pitch about the 

elastic axis. This simplified structural model, illustrated in Figure 2.1, serves as a demonstration of the technique and is described 

in many references (see Fung (1955)) as a common example. 

the equations are developed about an arbitrary static angle of attack, 

9, which is included in the aerodynamic model. 

However, 

The wing is attached to the support by linear springs to model 

stiffness characteristics and by linear dashpots to model damping 

characteristics. The point of attachment is referred to as the elastic axis. The center of gravity is offset from the elastic axis by 

the parameter r. This offset affects sensitivity to flutter; this 

effect has been described by Pines (1958) and Bisplinghoff (19751, as 

well as others. The aerodynamic model accounts for spanwise effects; 

hence, the aspect ratio of the wing is considered. Physical properties of the wing, stiffness (and structural damping) properties of 

the support, and aerodynamic properties act on the wing section as 

illustrated. 

17 
Figure 2.1. The Rigid-Wing Model. 
19 

2.2 Development of the Equations 

The components of the acceleration of the center of gravity are 

described in the inertial frame. The x-component is expressed as 

and the y-component is expressed as 

where a and y are the coordinates which describe the pitch and 

punge motion. The location of the center of gravity relative to the 

elastic axis is given by 

center of gravity is aft (i.e., towards the trailing edge) of the 

elastic axis. 

r, which is a positive quantity if the 

Summing forces in the X direction results in the following 

equation : 

Fx + N sin(a + 4) 

(2.3) 02 

= mE-& sin(a + 4) - a r cos(a + $11 

Summing forces in the Y direction results in the following equation: 

Fy - N cos(a + 4) 

*2 

= mi? + Er cos(a + 4) - a r sin(a + +)I 
20 

Summing moments about the wing's center of gravity results in the 

following equation: 

Mr + Ms + NG - rF cos(a + 4) Y 

- 

where x is the distance between the reference position of the aerodynamic loads and the mass center of the wing and is a positive 

quantity if the mass center is aft of the reference position. The 

reactions of the springs and dashpots are denoted by Fy, and 

Ms. 

and the result is the following differential equation: 

F,, 

Equations (2.3) and (2.4) are substituted into Equation (2.51, 

(2.6) .. 2 '0 I = NE + Mr + Ms - rN - mry cos(a + +) - r ma 

cg 

Equations (2.4) and (2.6) represent the differential equations 

which describe the motion of the system. Linear elastic spring and 

damping forces and moments are defined as follows: 

F = -c$ - kyy Y 

aa 

Ms = -C - 

a 

These expressions are substituted into Equations (2.4) and (2.6). 

Structural nonlinearities could be introduced in the spring and 

damping model. 
21 

The force and moment equations are now rewritten as 

.. 02 

my + mri cos(a + $1 - ma r sin(a + $1 

= -N cos(a + 0) - kq - Cy; 

2 .* I i + r ma + rmi cos(a + 4) cg 

. 2 .. = NX' - c a - kaa + Mr - rN - r ma 

a 

The form of these equations can be simplified by introducing the 

following: 

2 I =I +rm 

e cg 

and 

M = Mr + (Ti - r)N 

Equations (2.7) and (2.8) are now rewritten as 

I + mrf cos(a + 4) + kaa + ca & = M 

e 

and 

.. 

my + mri cos(a + $1 - mrG2 sin(a + 0) 

+ k y + cy; = -N cos(a + 0) Y 

(2.7) 

(2.10) 

The linear form of Equations (2.9) and (2.10) agrees with the 

equations derived by Fung (1955). 
22 

The numerical integration scheme which is described in Chapter V 

requires first-order differential equations. To this end, Equations (2.9) and (2.10) are solved simultaneously to provide equations 

in terms 

.. 

Y' 

.. 

a= 

of the second-order time derivatives. These equations are: 

(kaa + ca&r cos(a + $) 

1 

2 [Ie - mr2 cos (a + $1 

e 

I 

my 

+ lea .2r sin(a + $1 - - (k y + cy;) 

NIe cos(a + g) 

m 

- Mr cos(a + $1 - 

1 

- kaa + 

2 

- mr2 cos (a + $) ['e 

(2.11) 

2.2 - r a m sin(a + 4) cos(a + $1 

(2.12) 2 + r cos(a + $)(k + c i) + Nr cos (a + 4)) YYY 

where the aerodynamic normal force, N, and pitch moment, M, are 

defined as 

12 N * pUmACn (2.13) 

(2.14) 

M is referenced to the elastic axis. The quantity nllc is the full 

chord of the wing. 
2.3 Nondimenaionalized Form of the Governing Equations 

The governing equations of the UVLM have been written in a nondimensional form through the introduction of the characteristic 

length, tc, characteristic velocity, U,, and the resulting characteristic time, (T = tc/U,). 

dimensional by dividing both sides of the equation by 

equation is now expressed as 

Equation (2.11) is also rendered non2 U,/tCe This 

1 

K - r cos (a + 4) 

.. 

Y- 2 2 [ (Kaa + Ca& - nCCm)r cos(a + 4) 

+ K{i2r sin(a + 4) - K Y - C i - CC cos(a + +)}I Y Y n (2.15) 

2 where K = Ie/mtc 

C = pt A/2m 

Hence, we also utilize the mass, m, to characterize the physical 

properties in nondimensional form. Other terms in Equation (2.15) 

have been rearranged and represent the necessary collection of terms 

for the nondimensional form. All quantities are nondimensional. 

is introduced as a nondimensional inertia and 

is introduced as a nondimensional aerodynamic constant. 

C 

Similarly, Equation (2.12) is nondimensionalized by dividing both 

22 sides of the equation by Um/tc. This equation is now expressed as 

{<c i + K Y>r cos(a + 4) .. 1 a= 

Y 2 K - r2 cos (a + 4) Y 

(2.16) 
24 

One should note that the velocity now appears only in the spring and 

damping terms. K, and K,, contain the square of the speed in the 

denominator. C, and Cy contain the speed in the denominator. 

Hence, the effective stiffness and damping of the system decreases as 

the freestream speed increases. 
CHAPTER I11 

EQUATIONS OF MOTION FOR THE ELASTIC WING 

3.1 Description of the System 

In this chapter the equations of motion for an elastic wing are 

developed by using an energy approach. This wing is illustrated in 

Figure 3.1. Small displacements about an equilibrium position are 

assumed. The angle of attack is introduced into the equations. The 

wing is allowed to bend and twist about an elastic axis not necessarily coincident with the inertial axis. Thus, there is inertial as 

well as aerodynamic coupling between wing bending and torsion. 

The wing is represented as a cantilevered beam, which can both 

bend and twist. The cross section of the beam is assumed to be rigid. 

The wing root is rigidly fixed. The physical properties may vary 

along the span of the wing. Spanwise variation in the aerodynamic 

loading is considered. The solution of the governing equations consists of static and dynamic contributions. The statically deformed 

shape due to steady aerodynamic loads and the distributed weight is 

determined. Franz, Krenz, and Kotschote (1984) discussed the importance of deformations due to static loads on aerospace vehicles, 

particularly wind-tunnel models. 

3.2 Development of the Equations 

We ignore the out-of-plane (that is, the spanwise) motion. The 

components of the velocity are expressed as 

25 
26 

Figure 3.1. The Elastic-Wing Model. 
27 

v z = &y,t) + xa(y> cos(a(y,t)> &,t> 

The kinetic energy of the system is expressed as 

L1 

T = J 7 [m(y)({; + xa(y)(cos a>:]* 

0 

+ {xa(y)(sin 0);) 2 ) + Icg(y)i2] dy 

The potential energy of the system is expressed as 

(3.3) 

= L1 2 7 {EI(Y)W"~ + GJ(y)a' - m(y)g(w + xa(y) sin a)) dy (3.4) 0 

The distributed aerodynamic loads are provided by the unsteady 

vortex-lattice method (UVLM). The work done by the aerodynamic loading is expressed as 

(3.5) 

This aerodynamic spanwise loading is dependent upon the motion, 

and the motion, in turn, is dependent upon the aerodynamic loads. 

The normal force and pitching moment are €unctions which should be 

expressed as 

M = [w, G, w', a, G, t; 0, AR, planform1 
28 

According to Hamilton's principle, 

tL 

1 1 

00 

(6T - 6V + 6W> dy dt = 0 

Substituting T, V, and W into Equation (3.61, one obtains the 

following differential equations: 

.. *2 

mw + mx cos(a)& - ma sin(a)a a 

+ (E1 w")" - mg = -N 

The linear form of these equations is 

.. .. I a + mx w - (GJ a')' - mgxa = ?.I (3.10) e a 

where all properties vary along the length. 

For the case where Equations (3.9) and (3.10) represent a cantilevered wing, the boundary conditions €or these equations are 
29 

3.3 Static and Dynamic Equations 

The solution of Equations (3.9) and (3.10) provides both the 

static and dynamic contributions. The developed equations will yield 

the static solution for cases of decaying motion. However, the solution may require unnecessary computations for static solutions only; 

therefore, w(y,t) and a(y,t) are redefined as 

(3.12) 

(3.13) 

where the subscripts s and d refer to the static and dynamic solutions. Now these equations become 

(3.14) 1. + ma .. ad + (EIwi)" + (EIw'')" - mg = -N rnwd 

and 

I + mx - (GJa;)' - (GJa;)' - mgxa = M (3.15) ed ad 

From Equations (3.14) and (3.151, the equations for the static 

deformations and the dynamical equations for small motion about the 

static deformations are 

(EIw'')" - mg = -N 

S (3.16) 

S -(GJai)' - mgxa = M (3.17) 

(3.18) .. 

mwd + d + (EIw;)" = -N 
= Md I ed + muid - (GJai)' (3.19) 

Meirovitch (1980) describes a method of solution to these equations by 

an expansion of the dependent variables. The variables w and a 

are represented by expansions in terms of the natural free vibration 

modes of the system. 

I 

and 

J 

a = c qj Oj j =I+l 

(3.20) 

(3.21) 

where the qi and qj are the generalized coefficients and the Oi 

and Oj are the modes chosen as the comparison functions, and where 

I = number of bending modes selected, 

J - I = number of torsion modes selected, and 

J = total number of modes selected to represent the solution. 

3.3.1 Equations governing the static contribution 

These expansions are substituted into the static Equations (3.16) 

and (3.17). The result becomes 

I 

(3.22) 

3 - 

- 1 (GJO;)' qj = M~ + mgx, 

j =I+1 

(3.23) 
31 

Multiplying the first equation by ar for r=1,2,...,1 and then 

integrating along the span, one obtains 

LI L 

(3.24) 

Multiplying the second equation by (0, for S=I+l,...,J and 

then integrating along the span, one obtains 

L L 

(3.25) 

In Equations (3.24) and (3.25), the first term is integrated by parts 

and the boundary conditions, Equation (3.11), are imposed. With rearrangement these equations are rewritten as 

IL L L 

and 

JL L L 1 1 G.J@'@: dy q = MSas dy + 1 mgx, as dy 

j=I+1 0 j jo 0 

(3.26) 

(3.27) 

The above expressions represent J equations, which in matrix 

form may be written as 

(3.28) 
32 

The stiffness matrix [K] is partitioned as 

where 

L 

K12 = K = 0 21 

L 

= I GJ4'4' dy K22 ij 

for i = 1,2,.**,I 

j = 1,2,...,1 

for i = I+l,...,J 

j = I+l,...,J 

{As) is the static aerodynamic loading vector and is defined as 

for i=1,2,...,1 

j=I+l, . . . , J 

{SI is the static loading (due to weight) vector and is defined 

as 

for i=1,2,...,1 

j=I+l, . . . , J 
33 

3.3.2 Equations governing the dynamic contribution 

We now concentrate on the dynamic equations. Substituting the 

series representation into the governing equations ((3.18) and (3.19)) 

for the dynamic solution yields the following equations: 

I J J 

,! mXa Fi4, + 1 I ;i 4 - 1 qj(GJ4')' j = Md 

i= 1 j=I+1 e j j j=I+l 

(3.29) 

(3.30) 

Multiplying the first equation by 4r for r = l,2,...,1 and 

integrating along the span we obtain the following expression: 

I L J L 

I L L 

+ 1 qi / EX#;@: dy -/ Ndar dy (3.31) 

i=1 0 0 

Multiplying the second equation by (Ps for s = R+l,...J and 

integrating along the span we obtain the following expression: 

I L J L 

J L L 

(3.32) 
34 

The third term in each of these expressions has been integrated 

by parts and the boundary conditions, Equation (3.111, are imposed. 

The result is a more simplified form of the equations. The above 

expressions represent a number of equations (J) equal to the total 

number of modes chosen to represent the system. These equations may 

be written in the following matrix form: 

The mass matrix [MI is partitioned as 

where the terms in [MI are defined as follows: 

L 

L =I I00 dy. M22 e i j 

(3.33) 

for i = 1,2,...,1 

j = 1,2,...,1 

for i = 1,2,...,1 

j = I+l,...,J 

for i = I+l,...,J 

j = I+l,...,J 

The stiffness matrix [K] is as previously defined. 
35 

The right-hand side contains the aerodynamic contribution and is 

expressed as follows: 

for i = 1,2,.**,1 

j = I+l,...,J 

Nd and Md represent the sectional normal force and moment referenced at the elastic axis and computed from the dynamic solutions. 

[MI and [Kl are symmetric. In addition, these equations may 

be inertially uncoupled if the sectional center of gravity and elastic 

axis are coincident such that x, = 0. 

still coupled through the aerodynamic loads. 

However, the equations are 

The equations are to be cast in a manner similar to that of the 

two-degree-of-f reedom formulation. Therefore, TI = is introduced 

which allows the governing equations to be written in state space 

form. 

The above expressions provide a set of first-order differential 

equations, the solution of which describes the motion of the elastic 

wing. The integration scheme of these equations will be described 

subsequentially. 
CHAPTER IV 

AERODYNAMIC LOADS, THE UNSTEADY VORTEX-LATTICE METHOD 

4.1 Overview 

The equations of motion and the associated integration scheme are 

formulated in a manner which permit the aerodynamic loads to be determined concurrently with the motion. In this chapter, a technique is 

described which models the aerodynamic loads using the unsteady vortex 

lattice method (UVLM). This method has been completely described by 

Konstadinopoulos (1981); hence, the technique as applied in this 

research will be briefly reviewed. 

The WLM is a numerical model of the three-dimensional flowfield, 

which can be used to treat arbitrary maneuvers of wings of arbitrary 

planforms, including highly swept delta wings which exhibit leading 

edge separation, multiple closely coupled lifting surfaces, and lowaspect-ratio planforms. It can also treat arbitrary angles of attack 

and camber as long as stafl or vortex bursting in the near wake does 

not occur. The method accounts for the nonlinear effects of the wakes 

adjoining the tips and trailing edge. 

The ability of the WLM to predict accurate spanwise distribution of loads has been demonstrated. In addition to the work of 

Konstadinopoulous, it is also worth mentioning other applications of 

the UVLM. The method was used to model aerodynamics of close-coupled 

lifting surfaces (Elzebda (1986)). Mook and Nayfeh (1985) demonstrated the method for high-angle-of-attack aerodynamics. Kobayakawa 

and Onuma (1985) applied the vortex-lattice method to model propeller 

36 
37 

aerodynamics. Thomas and Nerney (1976) implemented the vortex-lattice 

method (coupled with slender body theory) to predict the aerodynamics 

of wing-body combinations. Kandil, Mook, and Nayfeh (1976) applied 

the UVLM to aircraft interference problems. 

The approach assumes that the flow is incompressible and inviscid 

and does not separate on the wing, but that separation occurs along 

the sharp edges where the Kutta Condition is imposed in a steady flow. 

At each time step this vorticity, which forms the wake, is convected 

at the local particle velocity; thus, the position of and the distribution of vorticity in the wake are part of the solution. The wake 

contains the history of the flow. Consequently, the velocity induced 

by it at the present time reflects the previous motion. 

In the aerodynamic model the vorticity in the flow is restricted 

to a thin region around the lifting surface and its wake. The flow 

outside this region is considered irrotational. The wing and wake are 

modeled as a sheet of vorticity. The wing portion, commonly referred 

to as the bound-vortex sheet, is specified, and as a result, a finite 

pressure jump exists across it. The wake portion, commonly referred 

to as the free-vortex sheet, is not specified but rather is force-free 

and is formed as part of the solution in such a way that no pressure 

jump exists across the wake. 

4.2 A Description of the Method 

4.2.1 The wing representation 

The vortex-lattice representation of a typical wing is illustrated in Figure 4.1. The wing is modeled by a lattice of discrete 
38 

Figure 4.1. The Unsteady Vortex-Lattice Method. 
39 

vortex filaments. Each element of area on the wing is surrounded by a 

closed loop of constant-circulation vortex segments. These segments 

join at the nodes, where the loops make sharp (typically 90") turns. 

A node is the intersection of the segments and is represented by a 

heavy dot in the figure. Control points are located at the center of 

each element; two typical control points are shown in the figure 

(denoted by small crosses). 

The individual discrete vortex segments connecting the nodes are 

typically members of two different loops (the exceptions are the segments along the edges). Consequently, the circulation around an 

individual segment is the algebraic sum of the circulations around the 

two loops to which it belongs. In Figure 4.1 the r's represent 

circulations around the individual segments, and the G's represent 

circulations around the closed loops of segments. For example, 

referring to the figure, rl is Gg - G4. 

4.2.2 Kinematic considerations 

Two coordinate systems are defined: the fixed inertial frame and 

a moving frame fixed to the wing and aligned with the wing root and 

elastic axis. This would be a deforming coordinate system for the 

elastic-wing model. 

erence frame. As described in the chapters addressing the equations 

of motion, the dynamic model will be posed in terms of a translating., 

but nonrotating, coordinate system. 

The problem is posed in terms of the moving refThe position of a point may be described as 

+++ 

p=R+r 
40 

The velocity of this point is given by 

+ + 

where v is the velocity vector relative to the moving frame and w 

is the angular velocity of the moving frame. The components of the 

velocity vectors are expressed in terms of the moving frame. The 

components of the angular velocity in the moving frame may be expressed in terms of the derivatives of the Euler angles. 

description of the coordinate transformations and Euler angles for 

general motion is presented in Konstadinopoulos's work. The transformations required for this research involve only the pitch angle. 

These transformations are described in the chapter which presents the 

development of the governing equation of motions. 

A complete 

4.2.3 The Biot-Savart Law 

The equation that serves as the foundation of the UVLM is the 

Biot-Savart Law (see Karamcheti (1980)) which gives the velocity 

induced at point P (see Fig. 4.1) by an individual vortex segment. 

(4.3) 

The graphical representation of this equation is illustrated in Figure 4.1. I' represents the circulation associated with a vortex 

filament described by n and h is the perpendicular distance from 

the filament to the point P. 

+ 
41 

Equation (4.3) is used to establish the aerodynamic influence 

matrix and the local velocity at a point as affected by all bound and 

free vortex filaments. The velocity field generated identically 

satisfies the continuity equation by this equation. The continuity 

equation for an incompressible fluid is 

++ v*v=o (4.4) 

4.2.4 The boundary conditions 

The flow must satisfy the following boundary conditions. The 

disturbance velocity must approach zero far away from the lifting surface and wake (satisfied identically through the Biot-Savart Law), and 

the relative velocity normal to the lifting surface must vanish on the 

lifting surface, 

+ (? - cis) n = 0 on the wing (4.5) 

+ 

where Vls is the velocity of the lifting surface and n is the 

normal to the lifting surface. In addition, for an inviscid fluid the 

Kelvin-Helmholtz theorem, which requires that all vorticity be transported with the fluid particles when the pressure is continuous, is 

used to obtain the position of the force-free wake. This theorem may 

be stated as 

Dr 

-= 0 Dt (4.6) 

where r is the circulation around any arbitrary closed loop. 
42 

The vortex sheet representing the wing is approximated by a lattice. Thus, the no-penetration condition, Equation (4.51, is enforced 

at only a finite number of points, the control points of the bound 

lattice. These control points are located at the centroids of the 

elements, and the normal vectors are obtained by forming the cross 

product of the diagonals in each element. 

This lattice serves as a computationally expedient imitation of 

the boundary layers on the lifting surface and the free-shear layers 

in the wake. The velocity field induced by the vorticity is computed 

according to the Biot-Savart law; consequently, the velocity field 

satisfies the continuity equation and decays far from the wing and its 

wake . 

4.2.5 The calculation of the circulation 

The circulations of the loops, the Gi in Figure 4.1, are obtained from the no-penetration condition 

-b -b -+ N 

= (VIS - VWli ni for I = 1,2,...,N 1 Aij Gj j =1 (4.7) 

where N is the total number of the elements representing the wing, 

and the Aij 

control point of the i- (receiving) element generated by a closed 

th loop of unit circulation around the j- (sending) element. At the 

i- control point (Vis - Vw)i is the velocity of the lifting surface 

minus the velocity induced by the wake elements. 

represent the normal component of the velocity at the 

th 

th 
43 

For a rigid wing, the problem is posed in a moving coordinate 

AiJ system attached to a reference configuration of the wing and the 

are constant. For an elastic wing the A must be computed at each 

time step. 

ij 

4.2.6 The unsteady Kutta condition 

The Kutta condition requires ACp to be zero along the wing tips 

and trailing edge. It is satisfied at each time step by convecting 

the vorticity along the sharp edges of the wing into the wake at the 

local particle velocity. Each node is displaced according to 

++ 

Ar = v At 

+ 

where A; is the displacement of the node, v is the local velocity 

relative to the moving coordinate system, and At is the time step. 

The local particle velocity is computed from 

+++ ++ v=V - VA -w x r 

+ 

where V is the absolute velocity calculated from the Biot-Savart 

Law. To maintain the wake in a force-free position, each node in the 

wake is also displaced according to Equations (4.8) and (4.9). 

4.3 The Impulsive Start of the UVLM 

Initially, the wing is represented by a lattice. There is no 

wake. Before motion begins all circulations are zero. The instant 

motion begins the circulations in the bound vortex lattice change. A 

discrete vortex line forms along the wing tips and trailing edge. The 
44 

existence of this starting vortex is dictated by the requirement of 

spatial conservation of circulation. This is the vorticity that is 

convected into the wake from the sharp edges of the wing at the particle velocity. 

The entire aeroelastic model has been nondimensionalized by defining a characteristic length, &,, and characteristic velocity, U,. 

From these definitions the nondimensional time step, AT, is defined 

as Ec/Uw. As much as possible, it is desirable to have uniform elements in the lattices. To achieve this end, we choose the characteristic length to be the chord of the rectangular elements. At each 

time step the free-vortex system is convected at the local particle 

velocity, Equation (4.8). As a result the streamwise lengths of the 

wake elements are approximately the same as the bound elements. One 

might also note that an increase in the number of chordwise elements 

is accompanied by a corresponding decrease in the actual time step. 

4.4 Graphical Representations of the Flowfield 

In Figure 4.2 a flowfield predicted by the general UVLM is illustrated. The wing-tip vortex system is clearly evident. It strongly 

influences the velocity at the control points along the wing tips. In 

Figure 4.3 the flowfield is shown in different views. 

the strength of the vorticity fields of both the bound and free vortex 

sheets is shown. 

In Figure 4.4 

4.5 The Aerodynamic Loading 

The aerodynamic loading is determined by calculating the pressure 

jump across each individual element in the bound (wing) lattice. The 
45 

Figure 4.2. Wireframe of Wing and Wake Lattice. 
Figure 4.3. Vortex-Lattice Representations. 
Figure 4.4. Wing and Wake Vorticity for Wings of Different 

Aspect Ratios. 
48 

pressure is calculated from Bernoulli's equation for unsteady flows. 

The pressure difference between the upper and lower surface of element 

i is given by 

(4.10) 

+ 

where AV is the discontinuity in the tangential component caused by 

the local vorticity. The average circulation around the discrete 

vortex segments on parallel edges of an element is considered to be 

the circulation around a sheet of vorticity parallel to and between 

the same edges. The vorticity is concentrated into cores to facilitate computation of the velocity field. The average circulation is 

related to the jump in the tangential component of velocity perpendicular to these edges, across the thickness of the lifting surface. 

The reader is again referred to Konstadinopoulos (1981). 

The pressure distribution is multiplied by the area of the 

element to provide the elemental force perpendicular to the surface. 

These forces are integrated along the chord and span to provide the 

aerodynamic forces and moments. 
CHAPTER V 

THE NUMERICAL SOLUTION SCHEME 

5.1 Overview 

The equations governing the motion of the structure have been developed. In addition, the model that determines the motion-dependent 

aerodynamic loads has been described. The interaction between these 

aerodynamic loads and wing motion presents an interesting challenge 

since the loads are necessary to predict the motion, yet the motion is 

necessary to predict the loads. Therefore, an interacting numerical 

integration scheme has been developed, which determines the motion, 

loads, and wing/wake geometry simultaneously. The method is based on 

Hamming’s Predictor-Corrector Method as described in Carnahan (1969). 

This integration technique coupled with the nondimensionalized 

form of all governing equations results in a reliable arrangement of 

the problem. The characteristic length, ‘lc, is the chordwise length 

of a lattice element. One time step is represented by the time necessary for a fluid particle to convect approximately one 

predictor-corrector method does not subdivide a time step of integraE,. The 

tion as do other techniques such as Runge-Kutta. Hence, the chordwise 

size of the convected elements remains nearly uniform for all 

calculations of the loads. 

5.2 The First-Order Differential Equations 

The development of the two-degree-of-freedom model for the rigid 

wing and the multiple degree-of-freedom model for the elastic wing 

49 
50 

both result in a set of first-order differential equations as developed in previous sections. The set of first-order equations which 

describe the two-degree-of-freedom motion (see Eqs. (2.15) and (2.16)) 

can be written as 

. - nCC,)r cos(a + 4) 1 

n= 2 [was + CaE 

K - r2 cos (a + 4) 

- ccn cos(a + (PI)] 2 + K{S r sin(a + 41) - K Y - 

Y 

1 

K - r cos (a + 4) 

{(Cyn + K Y)r cos(a + 4) 22 Y 

i= 

2 + c(n cm + r cos (a + 4)~~)) 

where 

. 

Yon 

The equations of motion for the elastic wing are repeated for 

convenience. 

(5.3) 

(5.4) 

(3.34) 
51 

In general, we have a system of first-order differential equations of the form 

fl(t, Y1, Y2, ..., Yn) dY1 -0 

dt 

f2(t, Y1, Y2, .*e, Yn) dY2 -= 

dt 

(5.5) . 

. 

n 

dY 

dt -P fJt, Y1, Y2, .**, Yn) 

where n is 4 for the two-degree-of-freedom problem and n is twice 

the number of modes for the elastic wing problem. 

5.3 The Basic Predictor-Corrector Method 

In both sets of equations the right-hand side represents the derivatives of the dependent variables. The general predictor-corrector 

scheme requires the values of the dependent variables at the current 

and three previous time steps. In addition, the derivatives at the 

current and two previous time steps are required. With this information the computation of the state variables for the next time step can 

be calculated. The approach to obtaining the information necessary to 

start the procedure will’be described later in this chapter. The 

integration of one time step will be described when all the necessary 

information is known. 
52 

The state variables at the end of the current time step are predicted by the following equation: 

j = 1,2,..*,n 

where the subscript i is the index for the time step and the subscript 0 signifies the initial prediction of the state variables. 

The time interval, At, is unity in our formulation due to the nondimensionalization. The local truncation error, ej,i, is estimated 

at the end of each time step and this truncation error is used to 

modify the predicted values to y*. 

112 e I +- 

yj*,i+l,o Yj,i+l,O 9 j,i (5.7) 

Then y* is used to obtain the derivatives (f*). Next, the 

corrector" equation is used to obtain the updated state variables. ?f 

- 1 

'j,i+l,M '8 ["j,i yj,i-2 

The corrector equation is used iteratively (M is the index for 

this iteration) until convergence on the state variables is achieved. 

Hence, j,i+l,M-l' are computed from y 7, i+ 1 , M- 1 ' The derivatives, 

the corrector equation uses the current information. In addition, it 

is important to note that the aerodynamic model is used each time the 

state variables are predicted or corrected. 
53 

The wake is convected from the position generated at the end of 

the last time step into the new position. This wake, new state variables, and wing geometry are used to compute the loads. Therefore, 

the aerodynamic loads are current. 

Once convergence of the state variables and loads has been 

achieved the local truncation error is estimated. 

1 P- 9 - 

j 9 i+l 121 ('j,i+l,M Yj9i+l,M-l e 

A final modification is made to the state variables. 

P -e yj 9 i+l 'j,i+l,M j,i+l 

(5.9) 

(5.10) 

5.4 Starting the Predictor-Corrector Method 

The predictor-corrector procedure requires the derivatives of the 

dependent variables from the current and three previous time steps. 

It also requires the derivatives of these variables for the current 

and two previous steps. 

dete'rmine these values as a starter for the predictor-corrector 

technique. We elect to use a Taylor series expansion to establish the 

starting values because this approach only requires the aerodynamic 

loads at integral time steps. 

For the ,th equation 

Carnahan suggests a Runge-Kutta scheme to 

(5.11) 
54 

where the right-hand side of Equation (5.11) is the right-hand side of 

Equations (5.1) through (5.41, or (3.34). Then, Equation (5.11) is 

expressed in difference forms between the 1st and 2nd time step. 

Y -Y 

At = fn (t, Y1, Y2, ...s Yn) "2 nl 

1 

Y -Y 

or 

y = f (t, Y1, Y2, ..., y 1 * At + y n9 n1 n n1 L I 1 

Now recall, 

2 

Y(t + At) = Y(t) + At Y'(t) +$- Y"(t) + h.o.t. 

(5.12) 

(5.13) 

(5.14) 

where the prime indicates differentiation and h.0.t. represents the 

higher-order terms. Therefore, 

2 

y = Y + At Y; +LY" + h.0.t. 

"3 "2 2 "2 

or, considering Equation (5.11), Equation (5.15) becomes 

2 f' At 

2! + h.0.t . "2 

Y =Y +Atf + 

"3 "2 "2 

(5.15) 

(5.16) 

for 2At and considering relation (5.14) 

ynl Approximating from 

+ h.0.t. 2! Y (t + 2At) = Y + 2At YA + 

"3 "1 1 

(5.17) 
55 

or 

Y = Y + 2At f + 2At 2 f: + h.0.t. 

"3 "1 "1 1 

(5.18) 

If we multiply Equation (5.16) by 4 and subtract Equation (5.18) we 

obtain 

At - 2f At + h.0.t. 

"3 "2 1 "1 n 

3Y = 4Y (5.19) 

Next we make the approximation f' = f' ; then Equation (5.19) can be 

"2 "1 

rewritten as 

Y 

Y =-Y 4 --+($fn2--f "1 ) At + h.0.t 

"1 3 "3 "2 

In a similar fashion an expression for Y is developed 

"4 

Y 

y =--y 4 -- n2 + (4 f - 2 f ) At + h.o.t. 

3 3 n3 3 n2 "4 "3 

(5.20) 

(5.21) 

Hence, the starting values for the general predictor-corrector 

procedure are available. 

The integration of the equations of motion by the predictorcorrector method is the same for both the rigid wing on the elastic 

foundation and the elastic wing on the rigid support. However, the 

complete integration process includes the interaction with the aerodynamic model; hence, differences may exist in the complete integration process. 

in the wing-fixed (not inertial) reference frame. Therefore, the 

The convection of vorticity from the wing is performed 
56 

basic difference in the integration process is whether or not the reference system is moving and the associated effect on the convection of 

the flowfield at the local velocity. Interestingly, the two-degreeof-freedom rigid wing presents a more complicated integration process 

than the multiple-degree-of-freedom elastic wing since the elastic 

wing is assumed to be cantilevered from a fixed support and the rigid 

wing is attached to a rotating support. The complete integration 

procedures will be described. 

5.5 Integration Method for the Rigid-Wing System 

5.5.1 First-order convection theory 

The integration technique for the equations developed for the 

rigid wing is illustrated by the flowchart in Figure 5.1. First, new 

state variables are predicted. These variables represent both static 

and dynamic contributions. Next, the wake is convected to the new 

force-free position from the position generated at the end of the last 

time step and the aerodynamic loads are computed. Then, the state 

variables are correct6d. It is important to note that the aerodynamic 

model is used each time the state variables are predicted or corrected. Convergence of both the state variables and the aerodynamic 

loads is required. 

convected from that generated (and converged) at the previous time 

If convergence is not achieved the flowfield is 

step, new aerodynamic loads are determined, and the state variables 

are again corrected. Convergence is checked. If convergence is 

achieved then the conditions for the beginning of the next time step 

are known. 
57 

. 

Conditions - Predict (y, 9, a, a) at t=i 

> 

Convergence 

(Y, PI a, u) 

(Cn, CmI 

I UVLM I 

- - , 

- Compute C, & C, 

Conditions 

at t = i+l 

u Correct (y, y, a, d) 

Recall 3 Flowfield 

at t=i 

Figure 5.1. Flowchart of the Integration Scheme Using LowerOrder Convection Theory, Rigid-Wing Formulation. 
58 

5.5.2 Second-order convection theory 

A higher-order convection strategy was examined (see Fig. 5.2). 

In this approach, a two stage corrector process was developed. 

stage one, the motion was corrected and the aerodynamic loads determined for a wake in a fixed position. That is, the position of the 

wake was not recomputed each time the state variables were determined. 

Stage one was repeated to convergence of both the loads and motion. 

In stage two, the wake was convected by an iterative process to a converged position with fixed state variables. Then, the field points 

for the wake were convected at a velocity based upon the average of 

the corrected variables and the variables from the previous time step. 

The new position of the wake is checked for convergence with the previous position. A new two-stage cycle begins: 

for stage one (motion and loads determined) and the converged new 

motion and loads were then used again for stage two (wake position). 

Motion, loads, and the wake were required to converge before the 

integration process continued to the next integration time step. The 

computational time was increased five-fold; however, the predicted 

results showed no discernible differences over the lower-order convergence scheme. 

In 

this new wake was used 

5.6 Solution of the Equations for the Elastic Wing 

5.6.1 Static solution 

The solution of the equations for the elastic wing is handled in 

a somewhat different manner. The static solution of the elastic wing 

is required prior to the start of the dynamic solution. The analysis 
Figure 5.2. Flowchart of the Integration Scheme 1Jsing SecondOrder Convection Theory, Rigid-Wing Formulation. 
60 

of the static problem begins with the specification of a nominal angle 

of attack. The corresponding distributed steady aerodynamic loads are 

computed. This load and the weight are then used to calculate the 

bending and torsional deflections. At this point the wing has a new 

shape, but the loads still correspond to the old shape. New loads, 

corresponding to the new shape, are obtained next; then the corresponding new shape is obtained. The procedure is repeated until 

either both the shape and loads converge or aeroelastic divergence 

(where the aerodynamic forces exceed the elastic restoring forces) 

occurs. The procedure is illustrated in Figure 5.3. The elements 

with dots at their centroids are part of the bound lattice (wing); the 

others are part of the free lattice (wake). 

5.6.2 Dynamic solution 

The integration technique which provides the dynamic solution is 

illustrated in Figure 5.4. Similarities do exist with the rigid wing 

formulation. The wake is convected to the new force-free position 

from the position generated at the end of the last time step. Unlike 

the rigid-wing formulation, this convection is only performed once. 

The new state variables are predicted for the first integration pass 

and are corrected for all subsequent passes. These state variables 

only represent the dynamic contribution. They are added to the static 

contribution and are used through the modal expansions to generate the 

new wing shape. The wake, new state variables, and wing geometry are 

used to compute the loads. Therefore, the aerodynamic loads are 

current. It is important to note that the aerodynamic model is used 
61 

WING PROPERTIES and 

FREESTREAM CONDITIONS 

GENERATE FLOWFIELD 

I RESULTS1 

Figure 5.3. Flowchart €or the Determination of the Elastic 

Wing Deformations Due to Static Loads. 
62 

CONDITIONS 

att-i 

I I 1 

CONVECT FLOW 

from CURRENT WING 

h WAKE POSITION 

1 PREDICT f :} I 

8 CORRECT f } 

1 

COMPUTE 

AERODYNAMIC 

LOADS (AI 

I 

AERODYNAMIC 

LOADS (A) I 

I CONDITIONS at t - i+l 

1 6 

Figure 5.4. The Integration Process for the Equations Caverning 

the Dynamic Solution of the Elastic Wing. 
63 

each time the state variables are predicted or corrected. The state 

variables are corrected, new wing geometry is computed, the aerodynamic loads are updated, and convergence of both the state variables 

and aerodynamic loads are checked. If convergence is not achieved the 

scheme returns to the corrector stage. Convergence is required of 

motion and loads at the end of each time step prior to advancing. 
CHAPTER VI 

APPLICATION OF THE MODEL TO SPECIFIC EXAMPLES 

6.1 Introduction 

In this chapter, the aeroelastic model is used with several examples to demonstrate the technique. These examples include cases for 

the rigid wing on an elastic support and the elastic wing on a fixed 

support. These examples illustrate the ability of the technique to 

predict the aeroelastic behavior and the unsteady aerodynamic loads of 

a wing. 

6.2 An Example of Aeroelastic Behavior for the Rigid Wing 

First, we consider an example that is similar to one discussed 

by Fung. 

about the elastic axis is placed in a uniform steady flow. A crossA large-aspect-ratio wing that can only plunge and pitch 

sectional view of this wing is shown in Figure 2.1. 

Fung's model for the aerodynamic loads is limited to twodimensional flow. Our model includes the effect of the wing tips and 

the wake; hence, these predicted results account for three-dimensional 

and unsteady characteristics. The essential physical properties €or 

this example are given in Table 6.1 (these are values used by Fung). 

No structural damping is present in this example; therefore, 

damping of the motion results strictly from aerodynamic effects. 

Later, we will demonstrate the additional effect that structural 

damping contributes to the wing response. In addition, the elastic 

axis and the mass axis are coincident; hence, the pitch and plunge 

degrees of freedom are not coupled inertially. The coupling which 

64 
65 

Table 6.1. Properties of the Rigid Wing Example 

Wing area 60 ft2 

Wing chord 60 ft 

Elastic axis/mass axis offset 0 

Elastic axis location (% of chord) 50% 

Mass of wing 

Mass moment of inertia 

Support translational spring stiffness 

Support torsional spring stiffness 

Translational natural frequency .880 rad/sec 

Torsional natural frequency I. 552 rad/sec 

Density of air 

2 lb sec 

269 ft 

150630 lb see2 ft 

208.5 lb/ft 

363020 f t lb/rad 

.002378 lb sec2/ft4 
66 

occurs between the degrees of freedom is a result of aerodynamic 

effects. 

Two different freestream speeds are considered for the example 

and the density of air at sea level is used. These conditions are 

chosen to bound the critical dynamic pressure which we define as the 

dynamic pressure at which the motion neither grows nor decays. This 

critical dynamic pressure results in wing flutter. 

Typically, for the two-degree-of-freedom model, a simulation is 

performed by establishing the flowfield for a wing which is fixed at a 

static angle of attack. Next, the wing is released from this fixed 

position and given an initial disturbance. The response of the wing 

is examined; a decay of motion in both degrees of freedom indicates 

stability and a growth in motion indicates instability. 

6.2.1 Subcritical response without structural damping 

In Figure 6.1 we show the response of the wing to a small initial 

disturbance which has the following form: 

Y(0) * 0 G(0) = 0.01 

i(0) = 0.02 

and 

a(0) = 0 

The freestream speed for this case is small and results in a dynamic pressure which is less than the critical dynamic pressure. The 

displacements and velocities associated with both degrees of freedom 

clearly decay with time. The pitch and plunge displacements reflect 
67 

-01 

0 

I 

I y -.01 

-.02 

- -03 

Cm 

-.12 L 

.02 

.01 

a0 

- .01 

- .02 

.I 

0 100 200 300 roo 500 

TIME STEP TIME STEP 

Figure 6.1. Response of the Rigid Wing at a Subcritical 

Dynamic Pressure. 
68 

the static and dynamic contribution; hence, in time the motion will 

decay and will eventually converge to the static (which includes the 

effects of both aerodynamic and mass induced loading) condition. 

In addition, as one would expect, the frequency of oscillation 

for each degree of freedom closely matches the uncoupled frequencies 

of oscillation associated with the respective elastic springs. The 

aerodynamic normal-force and pitching-moment coefficients are also 

shown. This aerodynamic loading decays to the steady aerodynamic 

loads. 

As a matter of interest, in nondimensional form the frequencies 

are dependent upon the freestream speed; therefore, a change in the 

freestream speed will change the nondimensional frequency of oscillation. However, for subcritical conditions the frequency in physical 

time would remain nearly constant. Also, referring to the nondimensional form of the governing equations of motion (see Chapter 11, 

Eqs. (2.15) and (2.16)) one finds that the stiffness of the support 

spring is reduced as the freestream speed is increased - in fact, the 

square of the speed appears in the denominator. The nondimensional 

structural damping, if present, is inversely proportional to the 

freestream speed. 

6.2.2 Supercritical response to a small initial disturbance 

Next, we show the response of the wing to a dynamic pressure that 

is larger than the critical dynamic pressure. The initial conditions 

are the same (see Eqs. (6.1)) as those used in the previous case. The 
69 

time histories of both degrees of freedom and their time derivatives 

are shown in Figure 6.2. 

The initial disturbance grows and, as we will show, can be expected to reach a limit cycle eventually. A gradual transition occurs 

very early in the motion and is noticed in the response associated 

with the plunge degree of freedom. The motion for the two degrees of 

freedom changes from the natural frequencies in these modes to a common frequency, which is nearly the frequency of flutter in the limit 

cycle. This common frequency lies between the two natural frequencies 

and is closer to the natural frequency in pitch than to the one in 

plunge - a characteristic of the flutter phenomenon which has been 

well documented. The time histories of the aerodynamic normal load 

and pitching moment coefficients are also shown. It must be emphasized that the aerodynamic loads and wing motion are dependent upon 

each other and a form of the solution for the loads or motion is not 

assumed. At each time step, the integration process iterates to a 

converged solution, predicting both the flowfield and the motion of 

the wing simultaneously. 

6.2.3 Supercritical response to a large initial disturbance 

In Figure 6.3, we show the response of the wing to a large 

initial disturbance which has the following form: 

Y(0) = 0 ;(o) = 1.0 

and 

a(0) = 0 i(0) = 0 
70 

0 c 

-.lo 

Y 

-.15 

- 20 

I II 

.2 r 

.I 

-.I 

- .2 t 

.02 .w 

.01 .02 

-.01 -.E 

- .02 -.w 

1 

. 04 

1 .o 

0 

-6 

Cn 

.2 

-.m 

- .2 

-.12 

1 1 1 1 1 

0 100 200 300 900 600 

1 1 IIIII 

0 100 200 300 yo0 500 

TIME STEP TIME STEP 

Figure 6.2. Response of the Rigid Wing at a Critical wnamic 

Pressure - Small Inftial Conditions. 
Y 

1.0r 

-1.oL 

-1.2 '*8B 0 100 200 300 9CO 5W 

TIME STEP 

.2 

a -2 

.lor 

.06 

.02 

a 

- .M 

- .06 

-.lOL 

T 

3 

1 

Cn 

Cm - .Y 01 -I 

-3 

100 200 300 900 500 

-5 

0 

TIME STEP 

Figure 6.3. Response of the Rigid Wing at a Critical Dynamic 

Pressure - Large Initial Conditions. 
The dynamic pressure is the same as the previous case. However, 

in this case we introduce a larger initial disturbance to only the 

plunge degree of freedom. 

The equations of motion are uncoupled in the absence of aerodynamic loads since r = 0. Yet, we find that the pitch is strongly 

excited almost immediately, which is clearly the result of coupling 

introduced by the aerodynamic loads. The plunging motion rapidly 

decays. In this case, the transition of the motion to a common frequency is clearly evident. The delay in this transition (when cow 

pared to the case with the small initial disturbance) is created by 

the large initial condition. 

After the initial rapid decay, the plunging motion appears to 

approach a limit cycle, not to decay to zero. The pitching motion 

continues to grow in time. Presumably, independent of the different 

initial disturbances introduced in these cases, the pitching and 

plunging motion of these two cases will eventually match. In Figure 6.4, we show the phase planes for the two cases shown in Figures 6.2 and 6.3. If we examine the phase planes associated with the 

plunge motion, we observe that in the upper half of Figure 6.4 the 

small motion appears to be growing around an unstable focus towards a 

limit cycle, while in the lower half of the figure the large motion 

appears to be decaying toward the same limit cycle. 

6.2.4 Subcritical response with structural damping 

No structural damping in the support is included in these previous results, though aerodynamic damping is modelled by the UVLM. 
73 

1 .o 

.2 

- .2 

E 1 

-1.01 

-5 -3 -1 I 3 5 7 

Y 

.10 r 

-.2 .2 

a 

Figure 6.4. Phase Plane Diagrams at the Critical wnamic 

Pressures - Small and Large Initial Conditions. 
74 

In fact, in some frequency domain solutions, the damping coefficient 

is used as an index to determine the onset of flutter. 

method can also account for damping in the support. In addition, the 

present equations of motion can be readily extended to include nonlinear stiffness and damping for both degrees of freedom thereby 

providing a means to study nonlinear structural behavior. 

The present 

In Figure 6.5, we show the response to the same initial conditions and same dynamic pressure as used in the case reflected in 

Figure 6.2. However, we now introduce structural damping into the 

problem, for the case illustrated Now instead of 

growing, the initial disturbances decay and the wing does not flutter. 

As one might expect, the flutter boundary is quite sensitive to structural damping. 

Cy = C, = 0.01. 

6.2.5 Comparison with theoretical results 

The unstable results discussed above were obtained by using a 

wind speed of 125 feet per second, and the stable results were obtained using a wind speed of 40 feet per second. The density of air 

at sea level was used and assumed constant. The critical speed lies 

between these two. Suggested methods for extracting the exact 

critical velocity will be discussed later. Fung found the critical 

velocity to be 162 feet per second which is higher than our predicted 

value. The differences in the two approaches are (1) Fung considered 

an infinite aspect ratio, and we consider an aspect ratio of ten which 

introduces finite wing effects; (2) Fung considered a zero static 

angle of attack, and we consider three degrees which introduce 
75 

.e 

.6 

Cn A 

.2 

0 

110 80 120 160 206 0 110 80 120 160 20; 

TIHE STEP TIRE STEP 

Figure 6.5. Response of the Rigid Wing at a Critical Dynamic 

Pressure With Structural Damping Present. 
76 

nonlinearities; and (3) Fung's aerodynamic analysis is based upon 

two-dimensional, linear theory, and our method is based on the 

unsteady vortex-lattice concept. 

6.3 An Example of Aeroelastic Behavior for the Elastic Wing 

We now turn our attention to results predicted for example cases 

of our elastic-wing model. These results describe the responses to 

initial disturbances for dynamic pressures below and above the flutter 

boundary . 

The physical properties of this elastic wing are given in 

Table 6.2. Although the structural model allows for the spanwise 

variation of all properties, we choose to demonstrate the procedure 

using a wing with constant sectional properties. The elastic axis 

is located at the sectional midchord, and the sectional center of 

mass is located aft of the elastic axis. The angle of attack is 

15 degrees, this angle is defined as the angle of attack prior to the 

introduction of loading due to mass, steady aerodynamic loads, and 

unsteady aerodynamic loads. 

The natural modes for the structure are required for the expansion of the dependent variables, w and a. This expansion is given 

in Chapter 111. The Hunter method, described by Gray (19871, is used 

to determine the natural modes for the elastic wing. This method 

solves the two-point boundary value problem by a transfer matrix 

approach to the finite difference equations. Of course, the modes may 

be derived by any available computational or analytical method. 
77 

Table 6.2. Properties of the Elastic Wing Example 

Wing aspect ratio 

Wing chord 

Bending stiffness 

Torsional stiffness 

Elastic axis/mass axis offset 

Distributed mass 

Distributed moment of inertia 

10. 

1.0 ft 

1.2 x lo5 lb-ft2 

7.0 x lo4 lb-ft2 

.1 ft 

.537 lb-sec2/ft2 

-125 lb-sec2 
78 

In this example we use the first and second bending modes and the 

a 

first torsion mode to represent the deformations of the elastic axis 

of the wing. The normalized amplitudes of the mode shapes are tabulated in Table 6.3. The natural frequencies, scaled to the fundamental frequency in bending, are l, 5.91, and 4.55. 

As discussed in Chapter V, the solution to the equations for the 

elastic wing consist of both the static and dynamic solutions. Tnitially, an undeformed wing is placed in a steady flow. The deformations due to the static loads and weight are computed. If conditions 

are such that the wing will aeroelastically diverge, this will be 

determined during the computation of the static deformations. The 

dynamic response will be computed about the statically deformed shape. 

6.3.1 Subcritical response of the elastic wing 

In Figure 6.6, we show the response of the elastic wing when the 

dynamic pressure is below the critical dynamic pressure. The initial 

disturbance has the following form: 

q1 = 2.5 

;Il = 0.0 . q2 f 0.0 

The generalized coordinates q1 and 92 are associated with the 

first two bending modes, 

mode. 

cients in the expansions for the dependent variables. 

43 is associated with the first torsion 

These generalized coordinates are the time dependent coeffiThese coefficients give the dynamic contribution only. 
79 

Table 6.3. 

Station 

(% wing) 

0. (root) 

01 

02 

.3 

04 

05 

.6 

.7 

.8 

.9 

1.0 (tip) 

Natural Modes Used to Represent the Elastic Wing 

Bending 

First Mode 

(Wbl = .044*) 

0. 

.0174 

.0650 

. 1380 

.2316 

.3413 

.4630 

.5926 

.7269 

.8632 

1 .o 

Bending 

Second Mode 

(wb2 = .260*) 

0. 

-. 0978 

-.3001 

-.5169 

-.6687 

-.6975 

-.5753 

-. 3064 

-0776 

.5287 

1.0 

Torsion 

First Mode 

(wtl = .198*) 

0. 

.1563 

.3088 

.4538 

.5875 

.7068 

-8088 

.8908 

.9509 

.9876 

1 .o 

*Frequencies in radians/nondimensional time. 
80 

Figure 6.6. Response of the Elastic Wing at a Subcritical 

Dynamic Pressure. 
81 

The motion decays for all three generalized coordinates to the 

statically deformed shape. In addition, the frequency of oscillation 

associated with the history of the generalized coordinates occurs at 

the natural, uncoupled frequency associated with each respective mode. 

6.3.2 Critical response of the elastic wing 

We now increase the dynamic pressure to that which creates the 

aeroelastic instability. This dynamic pressure was estimated by a 

method which will be discussed later. The initial disturbance has the 

following form: 

q1 = .25 q2 = 0.0 q3 = 0.0 

0 0 0 

q1 = 0.0 q2 = 0.0 q3 = 0.0 

These initial conditions are smaller than those used for the 

subcritical case. The motion predicted by the simulations for the 

critical case, which use the initial disturbance expressed in Equation (6.3), grew excessively fast due to the large initial deformations. Hence, smaller disturbances are used and reflect the same 

instability but at a lower initial rate of growth. 

In Figure 6.7 the response of the wing is represented by the time 

histories of the generalized coordinates. These coordinates are defined as in the previous case. The motion associated with the first 

bending mode (ql) decays initially; however, after this initial period 

the motion appears to neither grow nor decay. Several harmonics (from 

the different modes) are embedded in this motion. The motion associated with the first torsion mode grows in time. The frequency of 

oscillation for this mode is near the fundamental torsional frequency. 
82 

9' 

10 

6 

qa 

2 

-2 

-1 

0 200yoo8008001O0O 

TIE STEP 

Figure 6.7. Response of the Elastic Wing at a Critical 

Dynamic Pressure. 
83 

The history of 42 reflects the coupling necessary for flutter. 

Initially, the motion decays at the frequency associated with the 

second bending mode. Then, a shift in frequency and other harmonics 

appear. Around time step 600 one sees that the response (for all 

three coefficients) grows in amplitude. Finally, the frequency of the 

second mode coefficient coalesces with the torsional frequency. A 

flutter mode is created. The motion oscillates about a non-zero mean 

close to (but not coincident with) the static deflection. This drift 

is a nonlinear phenomenon (see Chapter 4 of the text by Nayfeh and 

Mook (1985)). 

6.3.3 Supercritical response of the elastic wing 

In the preceding case we considered a dynamic pressure very 

slightly above the critical dynamic pressure. We now increase the 

dynamic pressure to a larger value. The initial conditions are the 

same as those described by Equation (6.4). 

The response of the wing to these conditions is shown in 

Figure 6.8. An interesting feature is contained in the predicted 

response of the 'torsion coefficient, 93. The torsional contribution 

continues to grow. In fact, the oscillation grows to angles that 

would certainly cause the flow to separate from the surface, violating 

the assumptions used in the UVLM. However, more important is the fact 

that the motion is continually growing. Most likely, this motion , 

would be catastrophic to the structure. 

As a footnote, if we were to continue to increase the dynamic 

pressure it is possible that we could return to an aeroelastically 
84 

2 

1 

0 

-I 

0 -2 

TIHE STEP 

2 

1 

0 

-I 

0 -2 

TIHE STEP 

Figure 6.8. Response of the Elastic Wing at a Supercritical 

Dynamic Pressure. 
85 

stable regime. This characteristic is illustrated in Figure 1.3. As 

the freestream speed is increased, one or more modes may cross into 

the unstable regime. These modes may return to the stable regime (as 

does the mode 3 root) and the structure is again stable. For obvious 

reasons, the lowest critical condition is of interest in aircraft or 

wind-tunnel model design. 

6.3.4 Comparison with empirical results 

Harris et al. (1963) described a method to approximate the flutter speed based upon empirical data. 

flutter speed index adjusted according to particular characteristics 

(e.g. aspect ratio, mass ratio, sweep angle, mass distribution, etc.) 

of the configuration under study. These adjustments are based upon 

parametric analyses of experimental and analytical studies. 

The method uses a baseline 

For the elastic wing, the results reflect a constant freestream 

velocity which keeps the reduced (nondimensional) frequency nearly 

constant for all cases. The density of the fluid is modified which, 

in turn, modifies the dynamic pressure. 

Using the approach described by Harris, the density required for 

the system to become unstable is approximately 0.026 slugs per cubic 

foot. In our calculations, a density of 0.0066 slugs per cubic foot 

was used for the stable (subcritical) case, a density of 0.035 slugs 

per cubic foot for the unstable (supercritical) case, and a density 

of 0.030 slugs per cubic foot for the marginally unstable (critical) 

case. There is very close agreement between the empirical and 

numerically simulated condition for marginal stability. 
86 

6.3.5 Unsteady aerodynamic loads 

The aerodynamic loads and the motion of the wing are computed 

simultaneously and interactively. The primary objective of the 

numerical simulation is to determine the response of the wing and, in 

particular, determine whether or not the wing is aeroelastically 

stable. However, as a consequence, the unsteady aerodynamic loads are 

computed and the history of the pressure distribution on the wing 

could be used to describe the physics of aeroelastic response. Further, the technique could be modified by removing the structural model 

and restricting the wing motion to specified manuevers. 

In Figure 6.9 we show the pressure distribution of the wing at 

three different times during the wing motion. The spanwise pressure 

distributions are shown for several chordwise stations. The values 

of the pressure differences are for the control points of the wing; 

the control points along the wing tips, leading edge, and trailing 

edge are not located directly on the edges of the wing. 

In the upper portion the pressure distribution is shown for the 

wing due to static aerodynamic loads and the resulting static deformations. The pressure does decay toward the wing tips; however, a 

slight rise in pressure is indicated near these edges as a result of 

wing-tip vortices. The effect of the second bending mode is clearly 

evident. 

spanwise station) reflects that distribution one would expect from 

classical theory; the peak pressures are at the leading edge and these 

pressures decay to a near-zero value at the trailing edge as dictated 

The pressure distribution along the chord (for a constant 

by the Kutta condition. The pressure distribution is also shown at 

I' . 2 
87 

Steady flow 

o,o-o-o-o-o-o-o-o~o, 

0. 0, LEADING EDGE 

I 4 0 ACP 

0 TRAILING EDGE 

After 1 bending cycle 

r (1.5 torsion cycles) 

LEADING EDGE 

1 ACP 

0 TRAILING EDGE 

After 2 bending cycles 

(3 torsion cycles) f P 

LEADING EDGE 

1 AcP 

0 TRAILING EDGE 

0 .5 1 .o 

Spanwise location (% span) 

Figure 6.9. The Pressure Distribution of the Elastic Wing 

at Three Different Times. 
88 

two times during the motion of the wing. Again, the effect of the 

second bending mode is clearly evident. The effect of the strong 

wing-tip vortex systems created by the rapid motion is particularly 

evident in the bottom part of the figure. 

In Figure 6.10 the (nondimensional) time history of the pressure 

distribution for individual elements is traced. We show the pressure 

difference across four elements: the wing-tip leading-edge, the wingtip trailing-edge, the wing-root leading-edge, and the wing-root 

trailing-edge elements. The wing-tip elements experience largeamplitude motions; hence, large fluctuations occur. In contrast, the 

wing-root elements do not experience large motion because of the 

fixed-wing boundary condition; hence, smaller fluctuations occur. The 

pressure differences for the two trailing-edge elements are near zero 

which substantiates the presence of the Kutta condition in our 

aerodynamic model. 

In Figure 6.11 the (nondimensional) time histories of the aerodynamic normal-force and pitching-moment coefficients are shown for 

the wing-root, mid-span, and wing-tip locations. These coefficients 

are computed from the pressure distributions associated with the 

column of lattice elements located at these spanwise stations. As one 

would expect, the coefficients oscillate at frequencies associated 

with the natural modes. The peaks do not occur at the same time at 

all stations. 
AcP 

F I\ fl 

,/ \/LEADING EDGE, WING TIP 

A LEADING E., ROOT 

--< A TRAILING E., ROOT 

Time step 

Figure 6.10. The Predicted Pressure Distribution for 

Four Selected Elements. 
90 

root 

Cm 

s 

TIME 

Cn 

TIME 

root 

Figure 6.11. Sectional Aerodynamic Coefficients Predicted 

by the Model. 
91 

6.4 The Identification of Flutter 

We have shown the response of the rigid wing and the elastic 

wing to conditions which result in both aeroelastic stability and 

instability. However, the technique which determines the precise dynamic pressure necessary for the instability has not been discussed. 

Although not within the scope of this research, ultimately we are interested in determining the lowest possible dynamic pressure at which 

initial disturbances will not decay. 

ascertained then the effects of angle of attack, aspect ratio, or 

other nonlinearities may be examined. 

Once this exact condition can be 

Currently, we increase the freestream conditions (dynamic pressure) from a known stable position, simulate the response of the wing 

to an initial disturbance, and examine the response of the wing. At 

some point it is clear that the motion is growing and that a flutter 

mode exists. However, we have only bounded the critical dynamic pressure. It is difficult to determine the precise dynamic pressure at 

which the wing becomes unstable. In some cases, the motion associated 

with a degree (or degrees) of freedom appears to grow while the motion 

for another degree (or degrees) of freedom appears to decay. In 

addition, several simulations are required. A method is needed to 

identify this critical point with a minimum of required simulations. 

The frequency domain solutions (Desmarais et al., 1978) determine 

this critical point by tracking the eigenvalues (and associated aerodynamic damping) of the system for increasing velocities. We briefly 

describe these methods in Chapter I, the reader is referred to Figure 1.3. Hassig (1971) demonstrated a similar approach by examining 
92 

the damping associated with the assumed harmonic motion. Time domain 

solutions (Dowell et al. , 1980) typically use logarithmic decrement 

approaches to determine the flutter boundary. However, these methods 

assume a form of the solution. 

A method to determine the exact flutter boundary as predicted by 

the numerical simulation has not been fully developed; however, two 

promising approaches have been initiated. The first approach is 

based on a parametric identification procedure which is described by 

Konstadinopoulos (1984) and Elzebda (1986). The second approach is 

based on an energy method. The approaches will be briefly described. 

6.4.1 Parametric identification procedures 

Adopting the approach described by Konstadinopoulos, one can 

model the nonlinear motion of the wing by assuming the normal-force 

and pitch-moment coefficients are nonlinear functions of a, Y, and 

&. An equation for each aerodynamic coefficient is formed which 

contains the higherorder terms. A quintic nonlinearity is suspected; 

therefore, the equations contain all 55 terms up to and including the 

fifth-order terms. 

. 

The aerodynamic normal-force coefficient is represented by the 

following equation: 

02 c = al< + a a + a i + a 'i2 + a5a2 + a6a + ... n 2 3 4 

02 02 2.2 + a Y aa + a55Ya a + h.0.t. 54 
93 

The aerodynamic pitching-moment coefficient is represented by the 

following equation: 

-2 2 -2 C = blf + b2a + bjG + b4Y + b a + b6u + ... m 5 

-2 -2 2.2 + b Y aa + b Ya a + h.0.t. 54 55 

We then fit the aerodynamic coefficient versus time curves - such 

as those curves shown in Figures 6.1, 6.2, and 6.3 - using a leastsquares technique to determine the coefficients ai and bi. By 

taking off one term at a time we are able to identify those terms 

which have a negligible contribution. The remaining terms provide a 

simple expression for the aerodynamic coefficients. 

The coefficients ai and bi will be unique for a particular 

dynamic pressure, static angle of attack, spring constant, etc. 

Konstadinopoulos classified the coefficients as either restoring terms 

or damping terms. By examining the value and, in particular, the sign 

change of the damping terms versus the critical parameter (the angle 

of attack), he was able to identify the angle of attack at which wingrock occurs. 

Our critical parameter is the dynamic pressure and we are €nterested in the minimum value at which initial disturbances do not decay. 

This parameter will be identified by extracting from the damping 

coefficient versus dynamic pressure curve the condition at which the 

sign for the damping changes. 
94 

Unfortunately, our results have not yielded the same success as 

the work of Konstadinopoulos. Several complications are introduced 

when the method is extended to the aeroelastic application. First, 

we are working with two aerodynamic coefficients. Second, the terms 

of the fifth-order equations (see Eqs. (6.5) and (6.6)) which are 

eliminated are a function of the dynamic pressure; hence, the form of 

the expression is not quite the same for all cases. Third, the use of 

the method for the elastic wing will certainly be complicated by the 

spanwise variation and phase relationships of the aerodynamic 

coefficients. 

Yet, the concept of examining the character of the damping of 

the system is consistent with the frequency and time domain methods 

addressed in the literature. Therefore, the method may yet be the key 

in future determinations of the exact critical condition. 

6.4.2 Identification using the system energy 

An alternate method is suggested to determine the critical condition necessary for an instability. This method determines the total 

work being done on the system by the loads during the simulation. 

method is based on the concept that an aeroelastic instability exists 

when more energy is being extracted from the freestream by the wing, 

then pumped back into it by the wing. 

present when the total energy level of the wing grows in time. 

The 

Therefore, an instability is 

In Chapter I11 we develop the equations of motion for the elastic 

wing. The kinetic energy of the wing is expressed by Equation (3.3) 
95 

and the potential energy of the wing is expressed by Equation (3.4). 

These equations are repeated here for convenience. 

(3.3) 

V = 1 L1 7 [EI(~)w"~ + GJ(y)a' 2 

0 

The total energy is the sum of the kinetic energy and the potential energy. The state variables which have previously been predicted 

through a simulation and describe the motion of the wing are substituted into the energy expressions. 

Expressions that give the total energy for the rigid wing are 

available from the formulation of the equations of motion described in 

Chapter XI. The method is illustrated through an example utilizing 

the rigid-wing model. 

Figure 6.12 describes the time history of the wing in motion for 

both a stable and an unstable case. The predicted motion of the wing 

(Y,?,ct,&) is used in the expression for the total energy. The stable 

case is characterized by an oscillating, yet decaying, level of energy 

and will reach a constant energy level associated with the steadystate (i.e. static) conditions. The unstable case is characterized by 

an oscillating, yet growing, level of energy. If the dynamic pressure 
96 

Subcritical 

r I Critical 

TIME 

Figure 6.12. The Computed Energy of the Rigid Wing for a 

Subcritical and Critical Case. 
97 

is slightly above the critical condition then the motion may be 

governed by system nonlinearities and the total energy level will 

fluctuate around a nonzero mean. 

The advantage of the method is that now, instead of examining the 

motion associated with each degree of freedom or fitting the aerodynamic coefficients with high-order polynomials, a single quantity is 

used. 

mode 1 s . 

The method applies to both the elastic wing and the rigid wing 

The disadvantage is that a transition from the stable to the 

unstable case is still difficult to ascertain. One possible method is 

to measure the decrement associated with the decaying oscillation and 

extract a damping term from the energy curve. Negative damping 

indicates the instability. 

6.5 Wing and Wake Graphics 

An advantage of the general technique is that the time domain 

solution coupled with the WLM provides the capability of graphically 

depicting wing and wake motion. Now, the physics of aeroelastic 

response can be graphically described. 

We have previously demonstrated the capability (Chapter IV) of 

generating graphical representations of the flowfield. We now address 

the generation of movies of the wing and wake motion as predicted by 

the general model. Recently, we (Strganac and Mook (1987), and 

Strganac, Mook, and Mitchum (1987)) presented movies of the wing and 

wake motion. Both stable and unstable behaviors of the wing were 

shown . 
98 

In Figure 6.13 we show the vorticity for the wing at several time 

steps. The frames are extracted from the movie. The view is from 

above the wing surface. Near the top of each frame is the leading 

edge of the wing. The brighter regions, typically those regions 

near the leading edge, indicate the more intense concentrations of 

vorticity. As the wing twists and bends in response to the initial 

condition, the vorticity on the wing changes. It is important to 

reiterate that, in turn, the motion is also a function of this 

vorticity through the unsteady aerodynamic loads. 

In Figure 6.14 we show the vorticity and shape of the wing and 

wake at several time steps (which correspond to Fig. 6.13). The 

history of the previous motion is embedded in the wake. The character 

of the vorticity is also indicative of the current motion of the 

wing. The motion of the wing can also be seen in the sequence. 

For purposes of illustration, the physical properties of the 

system are such that for every three cycles of (first mode) bending 

motion there are two.cycles of (first mode) torsion motion. The 

geometry of the wing and wake continues to change. It is interesting 

to note that the wing-tip vortex system breaks down and reforms with 

the bending/torsion motion. The aerodynamic loads which are used to 

calculate the motion reflect this behavior. 
20th Time Step 

30th Time Step 

.( 40th Time Step 

Figure 6.13. The Strength of Vorticity for the Wing at 

Several Time Steps. 
Figure 6.14. The Predicted Shape and Strength of Vorticity 

for the Wing and Wake at Several Time Steps. 
CHAPTER VI1 

CONCLUDING REMARKS 

A numerical model has been described that predicts the steady and 

unsteady aeroelastic behavior of a wing. The method can predict 

steady-state static responses and transient responses to initial disturbances. The solution of the governing equations is obtained in the 

time domain. A physical interpretation of the response is aided by 

animation of the wing and flowfield histories. 

An aerodynamic model is introduced into the formulation which 

addresses nonlinearities created by static deformations, angles of 

attack, vorticity dominated flows, and unsteady flowfields. The wake 

is computed as part of the solution; hence, the history of the motion 

is included in the technique. Steady and unsteady aerodynamic loads 

are also computed as part of the solution. 

The numerical model has been developed such that the aerodynamic 

and structural models may be independently modified. As a result, 

nonlinear structural models or other aerodynamic models may be 

introduced in the future. Two structural models are demonstrated: 

the rigid wing mounted to a linear elastic support and the elastic 

wing cantilevered from a rigid support. Both formulations treat the 

nonlinear effects of static angle of attack and displacements. 

The equations governing the motion of the structure are coupled 

with the equations governing the motion of the fluid, and the structure and fluid are treated as a single dynamical system. The integration of the governing equations using the predictor-corrector technique requires convergence of the loads and motion at each time step 

of integration. 

101 
REFERENCES 

Ashley, H., Lehman, L. L., and Nathman, J. K., The Constructive Uses 

of Aeroelasticity, AIAA Paper No. 80-0877, AIAA International 

Annual Meeting and Technical Display, 1980. 

Atta, E. H., Unsteady Flow Over Arbitrary Wing Planforms Including Tip 

Separation, M.S. Thesis, Department of Engineering Science and 

Mechanics, Virginia Polytechnic Institute and State University, 

Polytechnic, 1976. 

Atta, E. H., Kandil, 0. A., Mook, D. T., and Nayfeh, A. H., Unsteady 

Flow Past Wings Having Sharp-Edge Separation, NASA SP-405, 1976. 

Atta, E. H., Kandil, 0. A., Mook, D. T., and Nayfeh, A. He, University 

Aerodynamic Loads on Arbitrary Wings Including Wing-Tip and 

Leading Edge Separation, AIAA Paper No. 77-156, 1977. 

Atta, E. H., Nonlinear Steady and Unsteady Aerodynamics of Wings and 

Wing-Body Combinations, Ph.D. Dissertation, Department of 

Engineering Science and Mechnics, Virginia Polytechnic Institute 

and State University, 1978. 

Batina, J. T., Seidel, D. A., Bland, S. R., and Bennett, R. M., 

Unsteady Transonic Flow Calculations for Realistic Aircraft 

Configurations, AIAA Paper Number 87-0850, 28th Structures, 

Structural Dynamics, and Materials Conference, April 6-8, 1987. 

Belotserkovskiy, S. M., Calculating the Effects of Gusts on an 

Arbitrary Thin Wing, Fluid Dynamics, Volume 1, No. 1, 1966. 

102 
103 

Belotserkovskiy, S. M., Calculation of the Flow around Wings of 

Arbitrary Planform in a Wide Range of Angles of Attack, NASA TT 

F-12,291, May 1968. 

Belotserkovskiy, S. M. and Nisht, M. I., Nonstationary Nonlinear 

Theory of a Thin Wing of Arbitrary Planform, Fluid Dynamics, 

Volume 9, No. 4, 1974. 

Bendiksen, 0. 0. and White, J. F. 111, Aeroelastic Tailoring of 

Advanced Composite Compressor Blades, AIAA Paper No. 86-1008, 

1986 . Biot, M. A. and Arnold, L., Low-Speed Flutter and Its Physical 

Interpretation, Journal of the Aeronautical Sciences, Volume 15, 

No. 4, April 1948. 

Bisplinghoff, R. L., Ashley, H., and Halfman, R. Lo, Aeroelasticity, 

Addison-Wesley Publishing Company, Inc., Reading, Massachusetts, 

1955. 

Bisplinghoff, R. L. and Ashley, H., Principles of Aeroelasticity, 

Dover Publications, Inc., New York, 1975. 

Carnahan, B., Luther, H. A., and Wilkes, J. O., Applied Numerical 

Methods, John Wiley d Sons, Inc., New York, 1969. 

Collar, A. R., The First Fifty Years of Aeroelasticity, Aerospace, 

February 1978. 

Crespo da Silva, M. R. M. and Hodges, D. H., Nonlinear Flexure and 

Torsion of Rotating Beams, with Application to Helicopter Rotor 

Blades - I. Formulation, and - 11. Response and Stability 

Results, Vertica, Volume 10, No. 2, 1986. 
104 

Desmarais, R. N. and Bennett, R. M., User's Guide for a Modular Flutter Analysis Softwear System (FAST Version loo), NASA TM 78720, 

1978. 

Devers, A. D., A General Method for Calculating Three Dimensional 

Nonstationary Aeroelastic Response in Subsonic Flows, Arnold 

Engineering Development Center TR-72-59, May 1972. 

Dowell, E. H. (Editor), et al., A Modern Course in Aeroelasticity, 

Sijtohoff d Noordhoff International Publishers, Alphen ann den 

Rijn, The Netherlands, 1980. 

Edwards, J. W., Applications of Potential Theory Computations to 

Transonic Aeroelasticity, ICAS-86 - 2-9.1, 1986. 

Elzebda, J. M., Mook, D. T., and Nayfeh, A. H., Unsteady Aerodynamic 

Interference for Lifting Surfaces, AIAA Paper Number 85-1801 , 

1985. 

Elzebda, J. M., Two-Degree-of-Freedom Wing Rock and Nonlinear Aerodynamic Interference, Ph.D. Dissertation, Department of Engineering 

Science and Mechanics, Virginia Polytechnic Institute and State 

University, December 1986. 

Farmer, M. G., A Two-Degree-of-Freedom Mount System with Low Damping 

for Testing Rigid Wings at Different Angles of Attack. NASA 

TM 83302, 1982. 

Franz, H. P., Krenz, G., and Kotschote, J., Load-Depending Deformations of Windtunnel Models, AIAA Paper No. 84-0589, AIM 13th 

Aerodynamic Testing Conference, March 5-7, 1984. 

Fung, Y. C., An Introduction to the Theory of Aeroelasticity, 

John Wiley & Sons, Inc., New York 1955. 
105 

Garrick, I. E., Aerodynamic Flutter, Volume V, AIM Selected Reprint 

Series, March 1969. 

Goland, M. and Luke, Y. L., A Study of the Bending-Torsion Aeroelastic 

Modes for Aircraft Wings, Journal of the Aeronautical Sciences, 

Volume 16, Number 7, July 1949. 

Gray, C. E. Jr., Coupled Bending - Torsion Steady-State Response of 

Pretwisted, Nonuniform Rotating Beams Using a Transfer-Matrix 

Method, M.S. Thesis, Old Dominion University, 1987. 

Hancock, G. J., Wright, J. R., and Simpson, A., On the Teaching of the 

Principles of Wing Flexure-Torsion Flutter, Aeronautical Journal, 

October 1985. 

Harris, G., Head, A. L. Jr., and Starr, G. A., Flutter Criteria for 

Preliminary Design, DDC No. AD-346962, 1963. 

Hassig, H. J., An Approximate True Damping Solution of the Flutter 

Equation by Determinant Iteration, AIAA Journal of Aircraft, 

Volume 8, No. 11, Nov. 1971. 

Hodges, D. H. and Dowell, E. H., Nonlinear Equations of Motion for the 

Elastic Bending and Torsion of Twisted Nonuniform Rotor Blades, 

NASA TN D-7818, 1974. 

Hsu, P. T., Some Recent Developments in the Flutter Analysis of 

Low-Aspect-Ratio Wings, Department of the Navy Contract No. 

NOa 55-771-c, 1957. 

James, R. M., On the Remarkable Accuracy of the Vortex Lattice Method, 

Computer Methods in Applied Mechanics and Engineering 1, 1972. 
106 

Kandil, 0. A., Mook, D. T., and Nayfeh, A. H., Nonlinear Prediction of 

the Aerodynamic Loads on Lifting Surfaces, AIAA Paper Number 

74-503, 1974. 

Kandil, 0. A., Prediction of the Steady Aerodynamic Loads on Lifting 

Surfaces Having Sharp-Edge Separation, Ph.D. Dissertation, 

Department of Engineering Science and Mechanics, Virginia 

Polytechnic Institute and State University, 1974. 

Kandil, 0. A., Mook, D. T., and Nayfeh, A. H., New Convergence Criteria for the Vortex-Lattice Models of the Leading-Edge Separation, 

NASA SP-405, 1976. 

Kandil, 0. A., Atta, E. H., and Nayfeh, A. H., Three Dimensional 

Steady and Unsteady Asymmetric Flow Past Wings of Arbitrary Planforms, AGARD Paper No. 2, AGARD Fluid Dynamics Panel Symposium, 

Unsteady Aerodynamics, Ottawa, Canada, 1977. 

Karamcheti, K., Principles of Ideal-Fluid Aerodynamics, Robert E. 

Krieger Publishing Company, Inc., Malabar, Florida, 1980. 

Kaza, K. R. V. and,Kvaternik, R. G., Nonlinear Aeroelastic Equations 

for Combined Flapwise Bending, Chordwise Bending, Torsfon, and 

Extension of Twisted Nonuniform Rotor Blades in Forward Flight, 

NASA TM-74059, August 1977. 

Kobayakawa, M. and Onuma, H., Propeller Aerodynamic Performance by 

Vortex-Lattice Method, Journal of Aircraft, Volume 22, No. 8, 

1985 . Konstadinopoulos, P., A Vortex-Lattice Method for General, Unsteady, 

Subsonic, Aerodynamics, M.S. Thesis, Virginia Polytechnic 

Institute and State University, July 1981. 
107 

Konstadinopoulos, P., Mook, D. T., and Nayfeh, A. H., A Numerical 

Method for General Unsteady Aerodynamics, AIM Paper Number 

81-1877, AIM Atmospheric Flight Mechanics Conference, August 

1981. 

Konstadinopoulos, P., Numerical Simulation of the Subsonic Wing-Rock 

Phenomenon, Ph.D. Dissertation, Department of Engineering Science 

and Mechanics, Virginia Polytechnic Institute and State 

University, March 1984. 

Konstadinopoulos, P., Thrasher, D. F., Mook, D. T., Nayfeh, A. H., and 

Watson, L., A Vortex-Lattice Method for General, Unsteady Aerodynamics, Journal of Aircraft, Volume 22, Number 1, January 1985. 

Konstadinopoulos, P., Mook, D. T., and Nayfeh, A. H., Subsonic Wing 

Rock of Slender Delta Wings, Journal of Aircraft, Volume 22, 

March 1985. 

Loring, S. J., General Approach to the Flutter Problem, Society of 

Automotive Engineers Journal, August 1941. 

Maddox, S. A., An Extension of Vortex-Lattice Method to Include the 

Effects of Leading-Edge Separation, M.S. Thesis, Virginia 

Polytechnic Institute and State University, 1973. 

Meirovitch, L., Computational Methods in Structural Dynamics, 

Sitjthoff & Noordhoff International Publishers, Alphen aan den 

Rijn, The Netherlands, 1980. 

Mook, D. T. and Nayfeh, A. H., Application of the Vortex-Lattice 

Method to High-Angle-of-Attack Subsonic Aerodynamics, SAE 

Technical Paper Series 851617, 1985. 
108 

Mook, D. T. and Maddox, S. A., An Extension of a Vortex-Lattice Method 

to Include the Effects of Leading-Edge Separation, Journal of 

Aircraft, Volume 11, No. 2, 1974. 

Nayfeh, A. H. and Mook, D. T., Nonlinear Oscillations, John Wiley d 

Sons, Inc., New York, 1979. 

Nayfeh, A. H., Mook, D. T., and Yen, A., The Aerodynamics of Small 

Harmonic Oscillations around Large Angles of Attack, AIAA Paper 

NO. 79-1520, 1979. 

Pines, S., An Elementary Explanation of the Flutter Mechanism, 

Proceedings of the National Specialists Meeting on Dynamics and 

Aeroelasticity, November 1958. 

Rivello, R. M., Theory and Analysis of Flight Structures, McGraw-Hill 

Book Company, New York, 1969. 

Strganac, T. W. and Mook, D. T., The Application of the Unsteady 

Vortex Lattice Method to the Nonlinear Two Degree-of-Freedom 

Aeroelastic Equations, AIAA Paper No. 86-0867, 27th Structures, 

Structural Dynamics, and Materials Conference, May 19-21, 1986. 

Strganac, T. W. and Mook, D. T., A New Method to Predict Unsteady 

Aeroelastic Behavior, AIAA Paper Number 87-0736, 28th Structures, 

Structural Dynamics, and Materials Conference, April 6-8, 1987. 

Strganac, T. W., Mook, D. T., and Mitchum, M. V., The Numerical 

Simulation of Subsonic Flutter, AIAA Paper Number 87-1428, 19th 

Fluid Dynamics Conference, June 1987. 

Theodorsen, T., General Theory of Aerodynamic Instability and the 

Mechanism of Flutter, NACA Report No. 496, 1940. 
109 

Thomas, J. L. and Nerney, B., Aerodynamics of Arbitrary Wing Body 

Combinations with Vortex Lattice and Slender Body Theory, AIAA 

Paper No. 76-198, AIAA 14th Aerospace Sciences Meeting, January 

26-28, 1976. 

Thrasher, D. F., Mook, D. T., Kandil, 0. A., and Nayfeh, A. H., 

Application of the Vortex-Lattice Concept to General, Unsteady 

Lifting - Surface Problems, AIAA Paper No. 77-1157, 1977. 

Thrasher, D. F., Nonlinear Unsteady Aerodynamics with Application to 

Dynamic-Aerodynamic Interaction, M.S. Thesis, Department of 

Engineering Science and Mechanics, Virginia Polytechnic Institute 

and State University, 1979. 

Weisshaar, T. A., Aeroelastic Tailoring - Creative Uses of Unusual 

Materials, AIAA Paper 87-0976-CP, 28th Structures, Structural 

Dynamics, and Materials Conference, 1987. 

Yates, E. C., Jr. and Bennett, R. M., Analysis of SupersonicHypersonic Flutter of Lifting Surfaces at Angle of Attack. 

AIAA Paper No. 71-327, 12th Structures, Structural Dynamics, 

and Materials Conference, April 19-21, 1971. 

Yates, E. C., Jr., Flutter and Unsteady-Lift Theory, Performance and 

Dynamics of Aerospace Vehicles, NASA SP-258, 1971. 

Yates, E. C., Jr., Wynne, E. C., and Farmer, M. G., Measured and 

Calculated Effects of Angle of Attack on the Transonic Flutter of 

a Supercritical Wing, AIAA Paper No. 82-0647, 23rd Structures, 

Structural Dynamics, and Materials Conference, May 10-12, 1982. 

Yates, E. C., Jr., Unsteady Subsonic and Supersonic Flows - Historical 

Review; State of the Art, Computational Methods in Potential 

Aerodynamics, Springer Verlag, 1985. 
VITA 

The author was born in Tucson, on He 

- 

attended high school at Frederick Military Academy in Portsmouth, 

Virginia and graduated in 1971. 

University from 1971 to 1977 and graduated with a Bachelor of Science 

degree in Aerospace Engineering. While at N. C. State University he 

participated in the Engineering Cooperative Education Program in which 

he was employed by the National Aeronautics and Space Administration 

He attended North Carolina State 

at Wallops Flight Center in Virginia. Upon his graduation he joined 

the staff of Wallops Flight Center on a full-time basis. The author 

attended Texas A6M University during 1980 where he completed his 

degree requirements for a Master of Science in Aerospace Engineering. 

In 1982 the author transferred from Wallops Flight Center to Langley 

Research Center in Hampton, Virginia to pursue his interest in 

aeroelasticity. The author attended Virginia Polytechnic Institute 

and State University during 1984 and 1985 where he completed his 

course work, initiated his dissertation research, and satisfied the 

residency requirements for this degree. He returned to Langley 

Research Center to complete his research. 

Kathy, and their two children, Christopher and Kasey, currently live 

in Hampton, Virginia. 

The author, his wife, 

110