CFDBASED AEROSERVOELASTIC PREDICTIONS
ON A BENCHMARK CONFIGURATION USING
THE TRANSPIRATION METHOD
By
COLE H. STEPHENS
Bachelor of Science
Oklahoma State University
Stillwater, Oklahoma
1995
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
July, 1998
OKLAHOMA STATE UNIVERSITY
CFDBASED AEROSERVOELASTIC PREDICTIONS
ON A BENCHMARK CONFIGURATION USING
THE TRANSPIRATION METHOD
Thesis Approved:
/
Dean of the Graduate College
II
ACKNOWLEDGEMENTS
1 would first like to take this opportunity to thank my wife, Angie, and my two
children, Hunter and Samuel, for their patience, support, and encouragement during this
endeavor. Often enduring long days and late nights, Angie deserves a special debt of
gratitude for always being there as a source of Christlike love and inspiration.
AdditionaHy, I would like to thank my parents, Boyd and Cheryle, for instilling within
me the values of diligence, hard work, and integrity. Special thanks also goes to my
grandfather, "Gangad", for his generous support of my education over the years.
I also extend my sincere appreciation to my major advisor, Dr. Andrew S. Arena,
Jr., for his continual support and guidance not only in research and academics, but also as
a personal friend. The high standards and expectations held by Dr. Arena are fueled by a
genuine enthusiasm towards both teaching and research. Special appreciation is also
extended to my committee members, Dr. P. M. Moretti, and Dr. G. E. Young for
comments and suggestions throughout this process. 1 also with to thank and acknowledge
Dr Kajal K. Gupta of the NASA Dryden Flight Research Center for his continued
SUppOlt and sponsorship.
Finally, my gratitude also goes out to Mr. R. C. Scott and Mr. M. R. Waszak of
the NASA Langley Research Center for the additional insight and data provi.ded on the
BACT windtunnel model.
III
Chapter
TABLE OF CONTENTS
Page
I. fNTRODUCTION 1
1. 1 Background , I
1.2 Problem Definition 3
1.3 Research Objective 4
2. LITERATURE REVIEW 5
2.1 Structural Dynamics 5
2.2 Unsteady Aerodynamics Solver 6
2.2.1 Transonic Small Disturbance & Full Potential Equations 6
2.2.2 Eul.er and NavierStokes Equations 7
2.3 Modeling Surface Deformations 8
2.3.1 BodyFitted Coordinate Systems 8
2.3.2 Dynamic Meshes 10
2.3.3 ReMeshing L3
2.3.4 Surface Transpiration........................................ . 14
2.4 Transpiration Concept 17
2.4.1 Origins of Transpiration ]8
2.4.2 Application to Current Research..... 19
2.5 Benchmark Models Program. . 27
3. METHODOLOGY AND PROCEDURE..... 31
3. I STARS Modules 31
3.1.1 SOLIDS Module 32
3.1.2 CFDModule 32
3.1.3 Aeroelastic and Aeroservoelastic Solver 34
3.2 Implementation of the BACT Model into STARS 35
3.2.1 BACT SOLIDS Model Development.. 36
3.2.1.1 Structural Mode Shape Definition 38
3.2.2 BACT CFD Model Development.. 43
3.2.2.1 BACT Geometry Specification in STARS 44
3.2.2.2 BACT Grid Specification in STARS 49
3.2.2.3 BACT Boundary Condition Specification in STARS 52
3.2.2.4 Effect of the Dissipation Parameters in STARS _ 53
3.2.2.5 SteadyState Solution Convergence Criteria........... 54
IV
3.2.3 BACT Uncertainty Estimation 55
3.2.4 BACT Aeroelastic/Aeroservoelastic Development 61
3.2.4.1 TimeStep Definition in STARS 65
3.2.4.2 Modal Identification Technique 70
3.2.4.3 System Identification Technique 72
3.2.4.4 Control Law Development.. 75
4. RESULTS 79
4.1 Steady Results Without Control Surface Deflections 79
4.2 Steady Results With Control Surface Deflections 84
4.2.1 Steady Solutions for Transpired and Actual Control Surface
Deflections 85
4.2.2 Steady Solutions for Transpired Control Surface Deflections
Compared With Experimental Data 90
4.3 Aeroelastic Results 96
4.3.1 Unsteady Data for the BACT and NACA 0012 Wings Tested at
Langley 96
4.3.2 System Identification Parameters and Effectiveness 97
4.3.3 Flutter Prediction Using STARS 108
4.4 Aeroservoelastic Results Ill
4.4.1 Control Law Development.. 112
4.4.2 Control Implementation into STARS 113
4.4.3 Flutter Suppression for the BACT Wing Using STARS 115
5. CONCLUSIONS AND RECOMMENDATIONS 124
5.1 Conclusions 124
5.2 Recommendations 126
BIBLIOGRAPHY 128
APPENDICES 131
APPENDIX A 132
AI STARSSOLIDS Data File (NOPAPA.DAT) 132
A2 STARSSOLIDS Generalized Mode 1 Displacement
Definition (Plunge) 138
A3 STARSSOLIDS Generalized Mode 2 Displacement
Definition (Pitch) 144
A4 STARSSOLIDS Generalized Mode 3 Displacement
Definition (Control Mode) 150
A5 Program to Write Noda' Displacement Data into STARSSOLIDS
Format 156
APPENDIX B 158
B1 STARSCFD Geometry Data File (BACTDAT) 158
B2 STARSCFD Background Mesh Data File (BACTBAC) 187
v
B3 STARSCFD Boundary Condition File (BACT BCO) 189
B4 STARSCFD Parameter Control File (RACT CGNU) 190
B5 STARSUnsteady Scalars File (RACTSCALARS) 191
B6 Portion of STARSUnsteady Arrays File (BACTARRAYS) 192
VI
•
LIST OF TABLES
Table Page
Table 31: Spreadsheet Layout for Manual Input ofa Structural Mode 39
Table 32: Actual and Suggested Number of Mesh Nodes for STARS Volume 51
Table 33: Experimental Measurements in Structural Parameters 57
Table 34: Experimental Measurements in Structural Parameters 57
Table 35: BACT Model Parameters in STARS 65
Table 41: Nominal and Actual Parameters for Mach 0.77, a=O°, 0=2° 90
Table 42: Nominal and Actual Parameters for Mach 0.77, a=O°, 0=5° 90
Table 43: Nominal and Actual Parameters for Mach 077, a=O°, 0=10°. ..91
Table 44: Nominal and Actual Parameters for Mach 0.82, a=O°, 0=2° 93
Table 45: Nominal and Actual Parameters for Mach 0.82, a=O°, 0:=5° . . 93
Table 46: Nominal and Actual Parameters for Mach 082, a=O°, 0=10°... . 93
Table 47: Model Orders and Associated RMS Errors at Various Mach Numbers 98
Table 48: Sensitivity to xcg for the NACA 0012 & BACT Wings in STARS .. 109
Table 49: Percent Change in CJfluller for a 09% Shift in xcg (8.028" => 8.10") 109
VII

LIST OF FIGURES
Figure Page
Figure 21: Physical Coordinate System 9
Figure 22: Computational Coordinate System 9
Figure 23: Reference Grid for Deforming Mesh Algorithm 11
Figure 24: Maximum Pitch Angle (a,=15°) Using a Deforming Mesh Algorithm 12
Figure 25: Minimum Pitch Angle (a=15°) Using a Deforming Mesh Algorithm 12
Figure 26: Slight Surface Element Rotation 15
Figure 27: Illustration of the Transpiration Concept 16
Figure 28: 2 x 1 Plate CFD Mesh 19
Figure 29: Actual 2x 1 Plate Deformation 19
Figure 210: Steady Pressure Contours on the Deformed 2x 1 Plate at Mach 0.95 20
Figure 211: AGARD 445.6 Wing, Undeflected and Deflected CFD Meshes 21
Figure 212 Steady Pressure Contours for the AGARD Wing at Mach 0.678 22
Figure 213: Variable WingFlap Intersection Example 23
Figure 214: Desired Flap Deflection......................................... 25
Figure 215: MeshShearing Example 26
Figure 216: Equivalent Mesh for Transpiration 26
Figure 217: BACT Wing ModeJ Dimensions 28
Figure 2J8: BACT Model on Flexible PAPA Mount.. 29
Figure 31: Stages of Advancing Front Technique 33
Figure 32 Block Diagram of TimeMarching Approach 35
Ylll

Figure 33: Finite Element Solids Mesh for BACT Wing 37
Figure 34: Rigid Body Pitch Mode Definition Example 40
Figure 35: Structural Deformation Opposite Original Definition .41
Figure 36: Implemented Rotational Mode Shape for STARS 42
Figure 37: CFD Computational Volume Specification 45
Figure 38: Wing Geometry Specification 46
Figure 39: CloseUp ofDeflected Flap Geometry Definition 47
Figure 310: Views Showing Tetrahedral Surface Mesh on the BACT Wing 50
Figure 311: Effect of Dissipation Constants on Cp in STARS CFD Solution 53
Figure 312: Solution Convergence Using Residuals and Maximum Mach Number. ..... 55
Figure 313: Effect ofKh, Ka, and xcg Location on Flutter Prediction 59
Figure 314: Effect of TimeStep at Two Different Values of ncycl on the Lift Evolution
for an Impulsively Started NACA 0012 Airfoil at Mach 0.3, a=5° 68
Figure 31S: Effect ofncycl with at Two Different TimeSteps on the Lift Evolution for
an Impulsively Started NACA 0012 Airfoil at Mach 0.3, a=5°. . 69
Figure 316 Abbreviated TimeHistory ofBACT Wing in STARS .. . 71
Figure 317: Complete SOOOStep TimeHistory ofBACT Wing in STARS 71
Figure 318: MultiStep Sequence for the BACT Wing (3Modes) 73
Figure 319: Modeled and Actual Response to MultiStep Input 74
Figure 320: MATLAB/SIMULINK® Model ofBACT with ControL 75
Figure 41: Steady Chordwise Pressure at Mach 0.51, a=O°, 0=0°, 60% Span .. g1
Figure 42: Steady Chordwise Pressure at Mach 0.67, a=O°, 0=0°, 60% Span. ... 81
Figure 43 Steady Chordwise Pressure at Mach 0 71 , a=O°, 0=0°, 60% Span .... ..... 82
Figure 44: Steady Chordwise Pressure at Mach 0.77, a=O°, 0=0°, 60% Span 82
Figure 45: Steady Chordwise Pressure at Mach 0.80, a=O°, 0=0°, 60% Span 83
Figure 46. Steady Chordwise Pressure at Mach 0.82, a=O°, 0=0°, 60% Span 83
IX

Figure 47: Comparison of Actual and Simulated 10° Control Surface Deflections 85
Figure 48: Surface Pressure Contours at Mach 0.77, 10° Control Surface Deflection ... 86
Figure 49: Comparison of Predicted Pressure Distributions for an Actual and
Simulated 10° Control Surface Deflection at Mach 0 77, 0° a 86
Figure 410: Pressure Contours at Mach 0.82, 10° Control Surface Deflection X7
Figure 411: Comparison ofPredicted Pressure Distributions for an Actual and
Simulated 10° Control Surface Deflection at Mach 0.82, 0° a 88
Figure 412: CrossFlow Velocity Vectors at the Trailing Edge ofthe BACT Wing
With an Actual and Si mulated 10° Control Surface Deflection.................... .... 89
Figure 413: Comparison of Predicted and Experimental Chordwise Pressure
Distributions at Mach 0.77, a=O°, 0=2° 91
Figure 414: Comparison of Predicted and Experimental Chordwise Pressure
Distributions at Mach 0.77, a=O°, 0=5° 92
Figure 415: Comparison ofPredicted and Experimental Chordwise Pressure
Distributions at Mach 0.77, a=O°, 0=10° 92
Figure 416: Comparison of Predicted and Experimental Chordwise Pressure
Distributions at Mach 0.82, a=O°, 0=2° ., 94
Figure 417: Comparison of Predicted and Experimental Chordwise Pressure
Distributions at Mach 0.82, a=O°, 0=5° SlS
Figure 418: Comparison of Predicted and Experimental Chordwise Pressure
Distributions at Mach 0.82, a=O°, 0= I0° . 95
Figure 419 Flutter Boundary Comparison Between 2 Geometrically Similar Wings:
NACA 0012 Wing (Air) & BACT Wing (R12) . 96
Figure 420: Training Data at Mach 0.51, q=14J.8 psf............................... . 100
Figure 421: Training Data at Mach 0.67, q=146.5 psf.. 101
Figure 422 Training Data at Mach 0.71, q=146.9 psL 102
Figure 423: Training Data at Mach 0.77, q=144.2 psf.. 103
Figure 424: Training Data at Mach 0.80, q=147.2 psf.. ., 104
Figure 425: Training Data at Mach 0.82, q=159.9 psf.. 105
x
Figure 426: Model Validation for STARS CFD/ASE: Model at Mach 0.&2 for Plunge
Displacement (in) _ 106
Figure 427: Model Validation for STARS CFD/ASE: Model at Mach 0.82 for Pitch
Displacement (deg) 106
Figure 428: Model Validation for STARS CFD/ASE: Model at Mach 0.82 for Plunge
Velocity (i nIsec) 107
Figure 429: Model Validation for STARS CFD/ASE: Model at Mach 0.82 for Pitch
Velocity (deglsec) 107
Figure 430: Flutter Boundary Prediction for Different xcg Locations 108
Figure 431 : STARS Flutter Prediction Compared with Experimental Results from the
NACA 0012 Wing and the BACT Wing _ 110
Figure 432: MATLAB® Flutter Suppression Example 112
Figure 433: Block Diagram of Control Implemented into STARS 113
Figure 434: Modified MultiStep Sequence Used for ASE 116
Figure 435: MultiStep Response for Model and Euler Solutions 117
Figure 436 Aeroservoelastic Response at Mach 0.51 119
Figure 437: Aeroservoelastic Response at Mach 0.77 120
Figure 438: MATLAB!l\l Model Comparison Using a Similar
Control Law at Mach 0.77 _.............. 120
Figure 439 Aeroservodastic Response at Mach 0.82 121
Figure 440: Euler Val.idation ofthe Modeled Aeroservoelastic
Response at Mach 0.77 _ 122
Xl

NOMENCLATURE
CFD Computational Fluid Dynamics
ASE Aeroservoelasticity
CG Center of Gravity
EA Elastic Axis
a Angle of Attack e)
8 Flap Deflection Angle (0)
M Mach Number
a Speed of Sound (in/s)
q Dynamic Pressure (psi, psf))
qf Dynamic Pressure at Flutter (psi, psf)
p Air Density (slinch/in3
)
y Ratio of Specific Heats (1.4 Air, 1.148 R12)
v Fluid Velocity (in/s)
0) Structural Natural Frequency (rad/s, Hz)
¢ Structural Mode Shape
COh Plunge Frequency (rad/s, Hz)
CO a Pitch Frequency (fad/s, Hz)
C08 Control Smface Frequency (rad/s, Hz)
Kit Plunge Stiffness (in Ib)
XII
m
10
Xcg
Pitch Stiffness (inolb/rad)
Mass (slinch)
Pitch Mass Moment ofInertia (slinchin2
)
Control Surface Mass Moment oflnertia (slinchoin2
)
Location of CG Relative to EA, Positive Aft (in)
PlungeP.itch Coupling (slinchojn)
PlungeControl Surface Coupling (slinchoin)
PitchControl Surface Coupling (slincho in2
)
Xlll
CHAPTER 1
INTRODUCTION
1. 1 Background
An efficient method of predicting the aeroservoelastic characteristics of modern
highspeed aircraft is crucial to aircraft design and flight testing. It is therefore essential
that the flight envelope be well defined prior to flight test operations. Without accurate
insight into an aircraft's aeroelastic tendencies, flight testing becomes a serious threat
both for the aircraft and its pilot.
Aeroelastic solutions are characterized by two maIO disciplines: structural
dynamics and computational fluid dynamics. Aeroservoelastic solutions include the
additional complexities introduced by forced control surface deflections during the
simulation. The structural dynamics portion of the code predicts a structures natural
response, or mode shapes. Any arbitrary deflection can therefore be described as a
superposition of a number of these natural mode shapes [Dowell, 1995]. Given an
arbitrary applied load, an aerodynamic load for example, the structural dynamics and
resulting deformations can be determined. The CFD solver uses the resulting
displacements and velocities that arise from the elastic structure and deflecting control
surfaces, and calculates new aerodynamic loads.

In the case of an aerodynamic body, these deflections have a great impact on the
flow field surrounding the body. Changes in this flow impact the lift, drag, and moment
experienced. This variation in loading is accompanied by a corresponding change in
structural deflections, which cause aerodynamic changes, which cause structural
deformations, and the cycle is repeated until one of two possibilities occur. One
possibility is that the changing aerodynamic loads and structural vibrations will
peacefully coexist and not result in a structural instability. The other possibility is that
the loads and deflections will coalesce and produce an unstable fluidstructure
interaction, also known as aerodynamic flutter. It is this flutter phenomenon that poses
the greatest threat to aircraft traveling at speeds ranging from high subsonic to
hypersonic. Allowed to progress, flutter has the definite possibility of causing structured
failure, and has the distinct probability of seriously injuring its pilot.
As described above, in the absence of forced control inputs, the classical
aeroelastic system simply reacts to the unsteady aerodynamics. In general, however,
aeroservoelastic systems have control surfaces such as ailerons and flaps that complicate
an aeroelastic analysis. Deflecting an aileron, for example, not only produces the
differential lift required to rolJ an aircraft, but also alters the twist of the wing itself. This
twist causes an effective increase or decrease, depending on how the aileron is deflected,
in the effective angle of attack seen by the entire wing. As a result, the effectiveness of a
deflected control surface decreases with increasing Mach number until the resulting
change in angle of attack exactly counteracts the increase or decrease in lift produced by
the aileron such that the aircraft does not roll. This aeroservoelastic phenomenon is
known as control surface reversal. In the case of flutter, control surfaces can serve as a
2
means by which to actively control aeroelastic response, falling under the category of
active flutter suppression.
Application of these solution techniques in an operational environment means that
the time it takes to complete a complete aeroelastic or aeroservoelastic simulation be kept
to a minimum wit.hout sacrificing solution accuracy. The structural solver requires far
less time, by several orders of magnitude, than does the CFD solution. Emphasis should
be given, therefore, to those means which improve the speed and efficiency of the CFO
solution.
1.2 Problem Definition
For current research, the STARS computer programs developed at NASA Dryden
Flight Research Center have been the primary means of a full ASE prediction [Gupta,
1997]. STARS is an highly integrated, finiteelement based code for multidisciplinary
analysis of flight vehicles induding static and structural dynamics, computational fluid
dynamics, heat transfer, and aeroservoelasticity capability.
Mentioned earlier, it is the CFD portion of the total simulation that requires the
vast majority of the solution time. Within each time step, structural deflections are
determined due to the predicted aerodynamic loads. Compared to the solution time
required by the CFD module, determination of the structural dynamics is essentially
instantaneous. This means that at each intermittent time step, it is the structural dynamics
solver that ends up waiting for the aerodynamic loads from the CFO portion of the code
This computational time is substantially increased if the solution must be paused at each
time step to deform the mesh based on a structural change due to modified aerodynamic
loads. Further difficulty is encountered if the mesh must be deformed in such a way as to
3
means by which to actively control aeroelastic response, falling under the category of
active flutter suppression.
Application of these solution techniques in an operational environment means that
the time it takes to complete a complete aeroelastic or aeroservoelastic simulation be kept
to a minimum without sacrificing solution accuracy. The structural solver requires far
less time, by several orders of magnitude, than does the CFD solution. Emphasis should
be given, therefore, to those means which improve the speed and efficiency of the CFD
solution.
1.2 Problem Definition
For current research, the STARS computer programs developed at NASA Dryden
Flight Research Center have been the primary means of a full ASE prediction [Gupta,
1997]. STARS is an highly integrated, finiteelement based code for mu.ltidisciplinary
analysis of flight vehicles including static and structural dynamics, computational fluid
dynamics, heat transfer, and aeroservoelasticity capabil ity.
Mentioned earlier, it is the CFD portion of the total simulation that requires the
vast majority of the solution time. Within each time step, structural deflections are
determined due to the predicted aerodynamic loads. Compared to the solution time
required by the CFD module, determination of the structural dynamics is essentially
instantaneous. This means that at each intermittent time step, it is the structural dynamics
solver that ends up wailing for the aerodynamic loads from the CFD portion of the code.
This computational time is substantially increased if the solution must be paused at each
time step to deform the mesh based on a structural change due to modified aerodynamic
loads. Further difficulty is encountered if the mesh must be deformed in such a way as to
3
account for discontinuous motions such as leading and trailing edge control surface
deflections. Accounting for these control surface deflections in a CPD grid presents
particular difficulty due the very close proximity of the control surface and adjacent wing
surfaces. In most cases, control surface deflections result in the exposure of surfaces not
previously seen by the CPD solver. These overlapping surfaces prove to be a significant
hindrance to flow computation.
] .3 Research Objective
In practical transonic and supersomc aeroservoelastic applications, thin,
Iightweight wings and control surfaces lend themselves to the susceptibility of flutter.
Along with continuing improvements in computational speed, there are more
sophisticated solution algorithms that take advantage of the additional speed and memory
capabilities. These advances in solution techniques continue to push the limits of even
the most powerful computers. In order to more fully appreciate advances in the stateof..
theart, ASE simulations must incorporate means which reduce the amount of
computational effort required to produce an accurate prediction. With the computational
overhead involved with timedependent deforming meshes, it is necessary to cultivate an
efficient means by which continuous surface deformations as well as control surface
deflections are accounted for in the ASE simulation as a means of actively controlling the
response of a system.
4
CHAPTER 2
LITERATURE REVIEW
Regardless of the solution methodology used, a fuIl ASE simulation requires a
means of coping with the structural dynamics and the determination of the natural mode
shapes, the unsteady aerodynamics, control inputs, and a means of incorporating these
structural and control surface deformations onto the CFD grid. Certainly, there will be
other differences within each simulation method, but at a minimum, the above items will
be common to virtually all ASE solutions.
2.1 Structural Dynamics
For the types of problems that are commonly encountered, the structural dynamics
portion of the solution is already much faster than the aerodynamics. The determination
of the structural mode shapes are generally determined one of two ways. First, the mode
shapes, pitch and plunge for example, could be known prior to the ASE simulation and
specified throughout the solution. A more general ASE simulation uses some sort of
structural dynamics solver, finite elements etc, to determine the structural characteristics
of the system. This type of solver computes arbitrary structural displacements based on
the aerodynamic loads. However, no matter how one chooses to solve for the structural
dynamics of the system, a significant amount of forethought must be given as to how
these structural deformations are related to a corresponding CFD grid. This point is
5


discussed in more detail later. STARS incorporates the finite element method to solve
for the structural response of the system.
2.2 Unsteady Aerodynamics Solver
The next issue is still the subject of a great many research papers. The question of
exactly how to model the unsteady aerodynamics is very often subject to computational
availability, time, and personal preference. Possibilities include, but are not limited to,
transonic small disturbance (TSD), and full potential equations (FPE), and more recently
Euler and NavierStokes equations. Historically, TSD and the full potential method were
most commonly used due to their compatibility with the computers of the time. With
advances in computer speed and memory, higher equation models such as the Euler and
NavierStokes have become more tractable.
2.2.1 Transonic Small Disturbance & Full Potential Equations
For three dimensional configurations, the transonic small disturbance equations
have been a popular choice for aeroelastic analysis and flutter predi.ction. The transonic
speed range is of primary interest because the flutter dynamic pressure is typically lower
there [Cunningham, Batina, & Bennett, 1988]. For the computational capability of the
day, the TSD and FPE equations were a popular choice because of their relatively low
computational cost and ease of implementation. Migration to more sophisticated models
is due mainly to the fact that these equations are not adequate in the presence of strong
shocks [Ruo & Sankar, 1987].
6
2.2.2 Euler and NavierStokes Equations
Advances in computational speed and memory have allowed the practical
implementation of Euler and NavierStokes solution algorithms to complex two and three
dimensional problems. These equations allow for the analysis of a wider variety of
problems at broader Mach number range than can be done with TSD or FPE equations.
The NavierStokes equations, with an adequate choice of turbulence model, are limited
only by the assumption that the fluid is a continuum. Take the viscous tenns out of the
NavierStokes equations, and the Euler equations are obtained. For sufficiently high
Reynolds numbers, the inviscid flow assumption makes good physical sense, as is shown
by the following equation:
Re = pUL = pU"
L f.J flU/L
(21)
The above equation expresses the Reynolds number as a ratio of the inertial forces to
viscous forces. It is apparent, therefore, that as the Reynolds number increases, ineliial
forces become more dominant than the viscous terms. The dominance of the inel1ial
terms in highspeed flows, such as those encountered during flutter, show that the
inviscid flow assumption made in the Euler equations are a valid means of aerodynamic
prediction. As one would expect, the Euler solutions are more limited in solutions where
there are significant boundary layer effects, boundaryIayerlshockintera.ctions, and
regions of separated flow.
Substa.ntial work has demonstrated the effectiveness of the Euler solution for
problems of practical interest. Free from the burden of determining a turbulence model
and constructing a mesh capable of resolving the boundary layer, an Euler solution is an
7
extremely attractive alternative to a code using the NavierStokes equations. Introduced
in section 1. 2, STARS makes use of the Euler equations on an unstmctured mesh for its
CFD prediction.
2.3 Modeling Surface Deformations
As with the choice of flow solvers, there are several popular methods of applying
a resulting surface deflection to a CFD grid. Many mesh deformation techniques use a
bodyfitted mesh which generally requires that the mesh move rigidly or shear as the
body deforms. These assumptions consequently limit the ASE analysis to rigidbody or
small amplitude motions [Batina, 1989]. Again, not an exhaustive collection of methods,
but a presentation of a few practical grid deformation techniques follows in the next few
sections.
2.3.1 BodyFitted Coordinate Systems
One popular method of accounting for stmctural deformations in the CFD mesh is
the use of a bodyfitted coordinate system. With this coordinate system, the wing surface
becomes a coordinate surface. This method involves a coordinate map from this physical
space to computational space [Malone, Sankar, & Sotomayer, 1984]. The relationship
between the physical and computational coordinate system can be visualized by
unwrapping the physical grid about a line, or axis, which lies within the wing surface.
Then, in the computational grid, the wing surface, as well as any assumed wake shape,
becomes a coordinate surface [Malone & Sankar, 1985]. Figure 21 shows the body
fitted coordinate system in the physical coordinate system. Note that key points are
8
labeled with letters. Figure 22 shows the transformed physical coordinate system in the
computationat coordinate system.
A
B
Figure 21: Physical Coordinate System
c
a ...  c' ~
Z \ Y\ P
x ~ TE \
N M
Figure 22: Computational Coordinate System
9


»
The above figures were obtained from a paper on the unsteady modeling of a fighter wing
in transonic flow [Malone, Sankar, & Sotomayer, 1984].
Resulting surface deformations must be related from physical to computational
space through a series of matrix transformations. These matrix transformations must be
calculated and implemented at each time step in an ASE solution. Although
implementation presents relatively few problems, the computational expense of these
transformations can be si.gnificant on complicated threedimensional geometries.
Additionally, this author has not seen this method implemented on a case involving a
discontinuous surface deflection such as those due to flaps or ailerons. As was discussed
earlier, the use of these meshes often require the assumption of smallamplitude, rigidbody
deformations.
2.3.2 Dynamic Meshes
Possibly the most intuitive of methods is the concept of a moving mesh. It simply
makes sense that one could deform the mesh in accordance to that predicted by a
structural dynamics solver. Work done by Batina has demonstrated the effectiveness of
such a method using an unstructured finitedifference mesh with an Euler solver [Batina,
1989].
Though the concept is simple, implementation comes at a price. What this type of
mesh boils down to is a large network of nodes connected by a series of springs whose
stiffness is inversely proportional to the length of its edge. At each time step, specified
boundary nodes are displaced by an amount corresponding to that of the aeroelastic
response of the body. The displacement of the rest of the computational domain is
therefore solved iteratively using static equilibrium equations in the x and y directions.
10

•
This results in x and y displacements for each of the interior nodes inside the
computational domain. This iterative procedure is accomplished by a predictorconector
method that first predicts the displacements due to linear extrapolation and corrects these
displacements with several Jacobi iterations of the static equilibrium equations.
Given in Figure 23, Figure 24, and Figure 25 are the original reference grid, the
deformed grid at maximum a and the deformed grid at minimum a, respectively, for a
wing oscillating about its quarter chord.
Figure 23. Reference Grid for Deforming Mesh Algorithm
Mentioned previously, the grid points on the outer boundary are fixed and the grid points
on the airfoil are fixed relative to the airfoil. From a maximum pitch oscillation of 15° to
a minimum pitch angle of 15°, the mesh smoothly transitions from one state to another
using the procedure described above.
11
rtf


Figure 24: Maximum Pitch Angle (a=15°) Using a Deforming Mesh Algorithm
Figure 25: Minimum Pitch Angle (a=15°) Using a Deforming Mesh Algorithm
The above figures were taken from a paper by J. 1. Batina [Batina, 1989].
12
,.
As one can imagine, the use of this type of mesh results in elements that have
been deformed from their original shape. These deformations lead to volumetric changes
within each element inside the computational domain. It is therefore necessary to add a
geometric conservation law to account for the changing cell areas at each time step. As
will be discussed later, deforming meshes also encounter difficulty in areas of surface
discontinuities.
Recently, an improved sprIng analogy was presented as an alternative to the
method proposed by Batina [Farhat, Degand, Koobus, and Lesoinne, ] 998]. In addition
to the linear springs between nodes, torsional springs at each node were also included to
further deal with the dilfficuhies involved with volumetric changes during mesh
deformation. Results were presented for a wing with a fulllength flap. Although related
to the problem of discontinuous surface deformations, the fulllength flap is more
amiable to this type of problem since moving surfaces never separate from one another.
Common to any dynamic mesh algorithm, substantial computational effort was involved
with deforming the mesh at each timestep. An esti mate was made that the

computational overhead involved in the implementation of this dynamic mesh accounted
for roughly 20% of the CPU time involved in a complete solution.
2.3.3 ReMeshing
Perhaps the most versatile option is the remeshing approach. Using this method,
the entire computational domain is remeshed at each timestep to account for structural
deformations and velocities. This method does not involve a complicated meshdeforming
algorithm, it simply redefines the surface geometry and generates a new
13


computational mesh. Of course, with current hardware, the remeshing approach is still
by far the most computati.onally expensive.
The problem with discontinuous surface deformations still exists with this
method. Even though the grid is redefined at each step and there is no meshshearing to
speak of, the varying intersection points at the interface of the wing and control surface
must still be calculated in order to model the geometry exactly. This calculation involves
specific knowledge about the geometry and would be difficult to implement in a generalpurpose
CFD code. Often, when a mesh is regenerated to account for control surface
deflections, additional surfaces are required to fill structural voids resulting from the
displacement. In an unsteady ASE simulation where the solution involves both wing and
control surface deflections, maintaining these varying intersection points would be
complicated at best and would most likely involve a substantial amount of user
intervention. Tills point is £luther illustrated in section 2.4.2.
23 4 Surface Transpiration
Though both the body fitted coordinate system and the dynamic mesh algorithms
have demonstrated their efficacy for solving aeroelastic problems, both require a
substantial amount of computational effort in between mesh deformations. As was seen
with the body fitted coordinate system, the resulting deformed grid must be mapped to a
computational system at each intermediate timestep. Even more so with the dynamic
mesh algorithms, consequential structural deformations result in a modification of the
entire computational domain.
In an environment where speed, without sacrificed accuracy, is of pnmary
concern, surface transpiration has shown itself as a viable tool to the aeroelastician. The
14

concept of surface transpiration is simple. With known stll.lctural displacements and
velocities, a simple modification to the nodal boundary conditions on the existing CFD
grid is capable of altering the displacements and velocities used in the flow solver.
With this method, no modifications are made to the existing CFD grid except for
a slight boundary condition modification to nodes on a deformable surface. As was
encountered with the previously discussed methods of grid modification, there are no
other complications associated with the transpiration method. With the transpiration
method there is no mapping from one coordinate system to another, no relative nodal
displacements, no elemental volume changes, no changes to the computational domain,
no need to iteratively solve for new nodal boundary conditions, etc. Stated again, the
only changes necessary are to nodal boundary conditions on deformable surfaces. Unlike
previous methods, deformations are accounted for only on those surfaces that require it.
What exactly is this change in existing boundary conditions? Generally it is quite
simply a change in the flow tangency boundary condition on an element. To attain the
noflow normal to the surface boundary condition, the flow solver computes a surface
normal for each surface element Observe Figure 26 below.
00r1>

Figure 26: Slight Surface Element Rotation
This figure shows an arbitrary surface element undergoing a sbght change in orientation.
It is important to keep the word slight in mind because it stands to reason that any
approximate method will loose effectiveness for large deformations. The figure shows a
15
single structural element with surface normal nl
being modified so its new surface
normal is n2 . Transpiration therefore assumes that there is no significant stretching or
volumetric change within the element so that the area of each element remains constant.
For a typical wing undergoing small amplitude structural deformations and control
surface deflections, this is a very reasonable assumption.
Assuming that a normal has an x, y, and z component, a change in orientation is
accomplished by changing the velocity boundary condition on the affected nodes. This
change in boundary condition comes in the form of an additional fluid velocity outside of
the existing surface elements. This additional velocity effects the way the unsteady flow
solver resolves the flow tangency boundary condition, see Figure 27 below:
Figure 27: Illustration of the Transpiration Concept
unsteady cases, the flow tangency boundary condition is represented by equation (22)
(23)
(22)
VTranspiral ion A
nOriginal
16
V·n =0
V·n =Vb·n
VOr;ginal
and (23), respectively.
the surface be deformed in such a way that it now has normal, nNew ' For the steady and
In the above figure, VOngmaliS the original tangential fluid velocity with normal, nOrlg;"<J"
Through an aeroelastic or control surface deformation, for example, the it is desired that


Equation (22) simply states that the velocity normal to the body must be zero Only
slightly more complicated, equation (23) states that the fluid velocity normal to the
surface must be equal to the velocity of the body nonnal to itself. In other words, no flow
can move through a solid surface. It is necessary to point out that the Vb mentioned here
is not the same as VTranspiration shown in Figure 27.
In summary, each surface element that is to undergo a change in orientation acts
as a source sheet. The strength of the source is determined by the extent of the simulated
deflection. Now, expand this procedure to an entire surface discretized into a large
number of elements. With a known surface deformation, perhaps from a finite element
solver e.g., it is desired that a surface be distorted from its original position. Within
reasonable limits, this arbitrary surface deformation can be simulated with an appropriate
change in the direction of the surface normal on each element making up the surface.
Since the flow solver is concerned with maintaining the flow tangency boundary
condition at each CFD node, the solution obtained on the simulated deformation should
closely approxi mate that of the actual deformation.
2.4 Transpiration Concept
Of the three methods of incorporating mesh modifications into the ASE solution
described in the previous section, the transpiration method shows the greatest potential
for accounting for mesh deformations with the least computational overhead. Its
simplicity is its greatest asset. Although the dynamics solver must still wail for the CFD
solver to predict the new aerodynamic loads, transferring the predicted deformations to
the CFD mesh is extremely fast. Since only the surfaces affected by the deflection are
affected, the rest of the computational domain remains untouched for the duration of the
17

ASE simulation. Surface normals on walls, farfields, and interior element surfaces are
also not modified. Appreciable time savings are realized due to the fact that a
modification to only those normals on the surface of a wing or fuselage, for example,
must be modified.
2.4. I Origins of Transpiration
Transpiration can trace its origins back to the late 1950's in a paper entitled On
Displacement Thickness which describes the "method of equivalent sources" for
modeling the influence of the boundary layer on the inviscid flow outside them [Lighthill,
1958]. Rather than thickening an actual airfoil, the boundary layer effect could be
accounted for by an equivalent surface distribution of sources. This is done by specifying
the necessary inflow or outflow boundary conditions on the original surface and solving
for the inviscid flow. As was described in Section 2.3.4, this method requires no
modification to the existing grid.
Simplicity, speed, and accuracy are the transpiration concepts greatest advantages.
As has been developed, the use of the transpiration boundary condition can be
implemented on an existing CFD grid with a minimal amount of computational cffOli.
The time it takes to simulate a deformed mesh is minimized due to the fact that no actual
grid deformation takes place, the computational volume is not modified, and only those
surfaces that require a boundary condition modification are affected. It's accuracy has
been effectively demonstrated over time through work done by Fisher, 1996, Raj &
Harris, 1993, Bharadvaj, 1990.
18
>'
2.4.2 Application to Current Research
Past research has demonstrated the effectiveness of the transpiration method when
applied to aeroelastic problems [Fisher and Arena, 1996]. For a vari,ety of problems
covering a wide range of Mach numbers, the transpiration method proved to be a viable
tool in the prediction of aerodastic responses. Here two specific examples are covered in
more detaiL The first is a 2x 1 plate case, the second is the AGARD wing.
The 2x 1 plate consists of a flexible plate surrounded by a rigid support, see Figure
28 below. To evaluate the usefulness of the transpiration method on this case, the CFD
mesh was deformed through a superposition of the first six natural modes, see Figure 29.
Figure 28: 2x 1 Plate CFD Mesh
Figure 29: Actual 2x 1 Plate Deformation
19
The transpiration method was used to simulate the actual deflection seen in the figure
above. For this case, at Mach 0.95, rdatively large surface deformations at this transonic
Mach number produce strong discontinuities on the pressures along the plate.
• e U8
·l.·.·~. 14 ~:~ \.r\J:" (j .{).2 : ' •
'{).4 ,
.{).6.:
.{).8~~1 'I ~.
·1 0 1 2 3
X
• ACtual
14 g:~ u~rOIl.'0J~'~
(j'()2j : I
'()4 J
.().6
.() 8 I~.'''l~_+~'+'''
~ 0 1 2 3
X
• c g~6. : rran.pifO""r\":'\... 02 :
14 0 •
(j '()2 . : •
.() 4 .• • •
'()6 .
.{)B ~~l~''~'<~'
~ 0 1 2 3
X
Figure 2 10: Steady Pressure Contours on the Defom1ed 2x 1 Plate at Mach 0.95
As can be seen in Figure 210, the transpiration method does an excellent job of
modeling the flow dynamics on the surface of the plate. In the figure above, three
lengthwise pressure cuts show the pressure distribution along each cut. In each section,
agreement between actual and simulated deflections are very good.
Another example of the application of the transpiration method is with the
AGARD 445.6 wing. This standard aeroelastic test case serves as a good reference for
application of the transpiration method to simulate surface deformations on a lifting
surface. Figure 211 shows two views of the AGARD wing. The leftmost figure shows
the undeformed mesh that will be used to simulate the figure on the right which is
20
actually deformed. This case serves to demonstrate the effectiveness of the transpiration
boundary condition when applied to relatively large surface deflection. As one can tell
from the figure, there are significant deformations resulting from both bending and
torsional modes.
Figure 211: AGARD 445.6 Wing, Undeflected and Deflected CFD Meshes
As was done with the 2x 1 Plate case, comparison is made between the simulated
and actually deformed mesh by means of chordwise pressure cuts at several points along
the span of the wing. For a Mach number of 0.678, we get Figure 212, below.
21

·1.2
1.0
0.8
0.6 8'OA
0.2
0.0
0.2
OA.
20
:H0j" .6· .• ATrcalnusapl iralJoo
_'''' ..... 8' 0.4 , L
.07" 0.2 :_==::::::::: 0.0 .• ~ ~
0.2. 
0.4 . ~"""I""""""""4~""""""""""'''''''.....o.....oL.L.''''''_~
10 15 20 25 30 35 40
X
:0.8~~G" • Aclual
0.6 . • Transpiration 8' 0.4 .'
02 . ~::=====::::=..... 0.0. • ::...
0.2· "'
0.4 '"t''''1......~...............~I...............
o 5 W U 20 ~ 30
X
Figure 212: Steady Pressure Contours for the AGARD Wing at Mach 0.678
For three chordwise pressure cuts through different spanwise locations along the wing we
once agam see excellent agreement between the simulated and actual surface
deformation.
What was lacking from the above two examples was a moving control surface.
Relatively smooth mesh deformations, as typically occur in aeroelastic problems, are
much more simple to deal with than are discontinuous surface deformations. For the
scope of the current research, the appealing characteristic about the transpiration method
is, oddly enough, the fact that the mesh does not move. Deflected control surfaces
provide several inherent difficulties for CFD solvers. When attempting to model a
control surface displacement, there are several factors that affect a CFD codes ability to
handle these difficult surface transitions.
22
First is the very close proximity of control surface edges to adjacent parts of the
airframe. Especially when using an Euler solver, these very narrow gaps present
significant computational difficulties. The flow through these gaps, along surfaces which
are parallel to the flow direction, will result in very high flow gradients and will
effectively wash out other, more significant, flow physics.
The second difficulty arises from the fact that even if one assumes that there is no
gap, the varying size of the face along the wingflap intersection would be terribly
difficult to account for, even in a dynamic mesh. Figure 213 below helps illustrate this
problem.
Figure 213: Variable WingFlap Intersection Example
Notice the area in the circled region in the above figure. For any change in flap angle, the
intersecting surfaces and the points of intersection change. Also observe that as the flap
changes position, the size and shape of the newly exposed surface changes. These
surfaces, specifically the lines defining the surfaces, must be modified with each different
£lap angle. The addition of these surfaces is necessary do keep the solution domain
closed. For the case of a wing with a finitespan flap, for example, deflection of the flap
requires the definition of 4 new surfaces with each new deflection. In either a dynamic
mesh or remeshing algorithm, for example, this variation in surface definition would be
difficult to account for.
23
, .1

Related to the second problem, is agam the difficulty encountered in the
immediate vicinity of the flap during a control surface deflection. With the flap in its
stowed position, there is essentially a smooth, continuous surface over the entire wing.
Assume that this flap, or control surface in general, is deployed several degrees. One
must consider what happens to the grid in the vicinity of the flap. With a dynamic grid,
the mesh must stretch to account for this displacement. The problem encountered with
this mesh deformation is the amount of mesh shearing that must be endured for the flap
to deflect.
Shown in Figure 215 is an example of this mesh shearing. For a simple wing
with a flap lying wilthin the span of the wing, a flap deflection similar to that of Figure 213
would produce surface discontinuities in the surrounding area of the flap. Figure 2] 4
shows the desired 10° flap deflection. The next figure, Figure 215, shows how a mesh
deforming algorithm might deform the existing mesh.
24
Figure 214: Desired Flap Deflection
Mesh shearing has the consequence of degrading the flow solution quality.
Notice that in the region of the flap, mesh shearing results in the elongation of elements
surrounding the wingflap intersection. Due to this shearing effect, there now exists poor
grid resolution around an area with high flow gradients, an area actually in need of grid
enhancement.
25
...... ;;;;;;;O;iiiiiiiiiiiiiiiill:l=_=:::::lO ;;;;......,;__._...
,r rELISAxPLT       §S~
Figure 215: MeshShearing Example
'), FELlSAXPLT B~ ~
Figure 216: Equivalent Mesh for Transpiration
26
_.' !.:J.
Using the concept of surface transpiration, there is no mesh deformation
necessary, hence no mesh to shear. Figure 216 shows the only mesh needed for the
application of a reasonably arbitrary flap deflection. With the above mesh, any arbitrary
flap deflection can be accounted for by simply rotating the elemental normals on the flap
by the desired flap deflection angle. Once again, one can see the speed at which this
method may be applied.
2.5 Benchmark Models Program
The Structures Division of NASA Langley Research Center (LaRC) initiated the
Benchmark Models Program (BMP) to obtain experimental data for the validation of
unsteady CFD codes. A variety of models were tested in the NASA Langley Transonic
Dynamics Tunnel (TDT) [Scott, Hoadley, Wieseman, & Durham, 1997]. In the BMP
program, two specific models are of interest. Each model has a rectangular planform
with a NACA 0012 crosssection, 16 inch chord, and 32 inch span. The first model was
simply a rigid rectangular wing fitted with pressure transducers over the sUliace of the
wing. The second model is referred to as the BACT, standing for Benchmark Active
Controls Technology. Though different models, each shares identical model dimensions,
and instrumentation. The only practical difference between the two models is the
presence of three control surfaces. These three control surfaces, two of which can be
seen in Figure 217, are a trailing edge control surface, and upper and lower spoilers.
27
...
);.
9.6"
.~"
./
I
../
"./
~.'
.1f
,I " I ,
19.2" ....
I I
"."
f
/
/
0'
,l2"
!
i'
f
I
i
10'

Figure 217: BACT Wing Model Di mensions
The control surfaces are centered along the models 60% span (19.2 in), and has a length
equal to 30% (9.6 in) of the wing's span. The trailing edge control surface has a width of
25% (4 in) model chord while the spoilers have a width of 15% (2.4 in) model chord.
The first model, the NACA 0012 wing, was tested in air and provided a large
experimental database. This database included steady pressure measurements, unsteady
pressure measurements during flutter, and flutter boundaries over a wide Mach number
range. Tested in R12, the BACT model's primary purpose was to provide additional
data for the purposes of evaluating a CFD code's effectiveness in modeling the control
surfaces illustrated above.
Both models were mounted inside the TDT on a device known as the Pitch and
Plunge Apparatus (PAPA) [Farmer, 1982]. Shown below in Figure 218, the BACT
model is seen mounted to the flexible PAPA mount system.
28
b
Figure 218: BACT Model on Flexible PAPA Mount
This mount system is simple and possesses dynamic properties that are easily obtained by
analytical means. It is important to note that the PAPA mount shown above is slightly
different than that described in the paper by Farmer, but these differences are primarily
cosmetic.
The mount basically consists of a model mounted to a "Chevron" bracket. Seen
on the Chevron mount are adjustable masses that allow adjustments to the models center
of gravity location. This Chevron mount is connected to a turn table by four steel rods
and a rectangular drag strut. The mount is designed such that it allows only two degrees
of freedom: rigid body pitch and plunge. The turntable allows an arbitrary choice in
angl e of attack. The Chevron mount, rods, drag strut, and turntable are hidden behind a
large splitter plate such that only the model is seen in the tunnel test section. For steady
pressure tests, this mount can be rigidified by replacing the Chevron mount, rods, and
drag strut by a large diameter (~6 in) rod.
With the quality and amount of experimental data available, these models serve as
the primary experimental benchmark to which all computational results obtained from the
29


current research are compared. Efforts presented within this paper illustrate the
implementation of the transpiration method within the STARS computer codes on:
steady pressure measurements, steady control surface deflections, conventional flutter,
and control inputs for purposes of flutter suppression.
30


CHAPTER 3
METHODOLOGY & PROCEDURE
The primary research tools for the current effort are the STARS codes developed
at NASA Dryden Flight Research Center [Gupta, 1997]. The current version of STARS
is the result of the evolution of the original STARS (STructural Analysis RoutineS)
computer code into an highlyintegrated multidisciplinary tool for the analysis of a wide
variety of 2D and 3D structures. This evolution involves the addition of several modules
to the original STARS code. Each individual module, general by design, is integrated
into an effective tool for the prediction of complicated aeroelastic and aeroservoelastic
problems. These modules include: structures, lI1eat transfer, linear aerodynamics, CFD,
controls engineering, and others.
3.1 STARS Modules
The scope of the current research is primarily involved with two of the modules
within the STARS computer programs. For a general ASE simulation, the user is
typically concerned with the structural dynamics of the system and the steady and
unsteady aerodynamic characteristics. The modules used for the current effort are the
structures and CFD modules, which are in turn integrated into the full ASE simulation.
31

3.1.1 SOLIDS Module
The SOLIDS module has a large solution bandwidth, but for problems pertinent
to current research, we are concerned with the determination of the free and forced
response. The fi"ee response comes from the solution of the foHowing equation:
3.1.2 CFD Module
where [M] and [K] are the inertial and stiffness matrices, respectively. Generally, once a
solids model is generated, STARS solves the above equation for the natural frequencies
(co) and mode shapes (4)). If, however, the natural frequencies and stmctural mode shapes
are known a priori, one can bypass this solution and manually create the generalized
mass and stiffness values.
The STARS flow solver is an Eulerbased code that applies finiteelement CFD
on an unstmctured grid. The implementation of an unstructured grid i.s a significant
feature of the STARS computer codes. For the general threedimensional case, the
computational mesh consists of an assemblage of tetrahedra. These tetrahedra are
oriented to form to the geometry being considered, thus making possible the treatment of
complicated shapes.
The unstructured grid shape is assembled using the advancing front technique.
This procedure consists of dividing a boundary into a finite number of points (nodes)
such that the external surface is sufficiently represented. Adapted from a figure by Peiro,
Peraire, and Morgan, Figure 31 shows how these triangles, or tetrahedra in three
dimensions, are arranged beginning at these outer nodes [Peiro, Peraire, and Morgan,
'1
[M]{ii} +[K]{u} =0 (31 )
32
. i

2
4
6
3
1
5
Figure 31: Stages of Advancing Front Technique
1993]. Additional tetrahedra are added in such a manner that the surface front collapses
upon itself until the entire domain is filled with tetrahedra.
The CFO module, in general, consists of four major parts:
• SURFACE: Generates the twodimensional front
• VOLUME: Generates the threedimensional computational domain
• SETBND: Defines the boundary conditions in the domain
33

order.
3.1.3 Aeroelastic and Aeroservoelastic Solver
[M]{u} + [c]{u }+ [K]{u }=f (t) (32)
[M] = generalized mass matrix
[K] = generalized stiffness matrix
[C] = generalized damping matrix
{u} = generalized displacement vector
• EULER: Steady or unsteady Euler flow solver
In general, the equations of motion for the coupled, timemarched ASE solution
dense arrangement of elements.
The user is able to specify certain parameters peliaining to the density of the CFD
adjust the existing computational grid such that regions of high gradients receive a more
involves the solution of (32), which is a matrix equation of motion for an arbitrary
Each one of the above steps, as one would surmise, need to be done in that particular
structure in general ized coordinates.
low mesh density in the farfield. STARS also has the capability of adaptive remeshing.
example, the user may wish to define regions of higher mesh density, while maintaining
surface and volumetric mesh. For regions such as leading and trailing edges ofwings, for
In the above equation:
Once a flow solution is obtained, the user has the option of letting STARS automatically
f(t) = generalized aerodynamic force vector
The general procedure, therefore, for solving aeroelastic and aeroservoelastic problems is
as follows. A steady CFD solution serves as the initial conditions for the structural
34
, J ._ ......
displacement and velocity boundary conditions serve as boundary conditions for the next
dynamics solver. A perturbation about this steady CFD flow will cause a change in the

structural displacement and velocity boundary conditions. These changes in
timestep in the CFD solution. Resulting forces and moments are then fed into the
structural dynamics solver which in turn computes new displacement and velocity
boundary conditions for the CFD flow solver. This process continues until the complete
timehistory is obtained. The above procedure can be visualized graphically througl1
Figure 32 below.
Global TimeStep
Figure 32: Block Diagram ofTimeMarching Approach
3.2 Implementation of the BACT Model into STARS
oos
FEM Solids SteadyState
Analysis CFD Solution
Modal ~ ! Initial Condili
Parameters 1 II Dynamics Solver I 1CFD Boundary Conditions
I CFD Solution I
I
Aero. Forces
As was mentioned in Chapter 2, the primary test cases for this effort are a NACA
0012 wing tested in air and the BACT wing tested in R12. Both in the Benchmark
Models Program and geometrically similar, models were tested under similar Mach
numbers and dynamic pressures in the Langley Transonic Dynamics Tunnel at NASA
35
Langley Research Center. The main difference, other than the fluid medium, is the fact
that the BACT wing has the capability of modeling control surface deflections, whereas
the NACA 00] 2 wing is simply a rigid, rectangular wing with no control surfaces.
The next few sections discuss in more detail, the incorporation of both these
models into STARS. Since, for all practical purposes, the two wings are exactly the same
except for the trailing edge control surface, the solids and CFD models used in STARS
will be the same..
3.2.1 BACT SOLIDS Model Development
Described in this section is an overview of the various steps taken to construct a
finite element solids mesh to represent the BACT wing. A solids model that included the
PAPA mount described in Chapter 2 was not developed due to the simple mode shapes
and frequencies exhibited by the BACTPAPA system. Since the model was constrained
to only plunge and pitch, the mode shapes and natural frequencies were available from
experimental data. Additionally, modeling the mount would have required a significant
amount of parameter finetuning in order to assure the natural frequencies and mode
shapes coincided with experimental results.
Even though the structural dynamics of the entire system were already known, a
structural mesh of the wing and flap itself is sti.1I needed. Mode shapes defined with this
model are in turn interpolated a CFD mesh. For too course a grid, it is possible that one
may introduce errors in the interpolation from one grid to another. For thi.s reason, care
was given to provide a tighter mesh in the region of the trailing edge control surface.
36
Figure 33 shows the resulting finite element structural mesh used in STARS.
rest of the wing.
f
Figure 33: Finite Element Solids Mesh for BACT Wing
rigid body pitch and plunge motions are encountered due to the PAPA mount. An
ICFD Case: NfA
Solids Case: bact
trailing edge control surface. To minimize any possible error in the interpolation from
obvious exception to the otherwise coarse mesh is the fight mesh in the region of the
Over the majority of the wing, a relatively coarse mesh is used due to the fact that only
the solids mesh to the CFD mesh this region was meshed much more densely than the

The majority of the solids mesh construction took place within the solids
preprocessor within STARS. PREPROCS is simply a tool that guides the user through
the creation of the mesh, assignment of structural properties, etc. It also formats and
37

writes the corresponding data file containing solution parameters, nodal properties and
locations, element properties and connectivity, structural properties, materials etc. The
thick lines running from leading to trailing edge along the edge of the trail ing edge
control surface and along the beginning of the flap actually correspond to small elements
that had to be added manually. It was discovered that in the interpolation from the solids
mesh to the CFD mesh exaggerated corresponding displacements of the flap due to the
large elements adjacent to the flap. Due to the way STARS implements the interpolation,
control surface deflections were seen out to locations corresponding to one half the size
of the larger elements. This was the first time thi.s problem had been encountered within
STARS. Before, when modeling continuous structural deformations, this sort of problem
never surfaced. To alleviate this problem, smaller elements had to be added around the
entire perimeter of the flap. Deflections still get interpolated out to one half of the
elements Width, but the elements are sized such that these errors are negligible.
In the PREPROCS routine discussed above, elements are assigned material types
and associated constants, nodes are constrained, etc. These constants, however, are not
used in this particular case for reasons discussed previously. The input data file created
by PREPRoes is given in Appendix A1. Since the structural mass, damping, and
stiffness characteristics are obtained experimentally, STARS allows the manual input of
this data. This point is covered in more detail a little later.
3.2.1.1 Structural Mode Shape Definition
In general, once one has developed the STARS solids mesh, the solution can be
set up to run the undamped, free vibration analysis to determine a user defined number
of structural natural frequencies and mode shapes. For the case of the BACT, the two
38
,
•

structural mode shapes are well defined, as is the control mode, and it is a relatively
straightforward procedure to define ones own structural mode shapes. Discussed in the
next few paragraphs is a general set of steps used in the creation of the two structural
mode shapes and the control mode.
First, the necessary parameters in the soEds file are set to solve for the first three
natural mode shapes. Run the undamped, free vibration analysts to obtain a properly
formatted out. 2 file. Although the data in the file will be replaced, STARS requires
proper formatting fOf later modules so the file serves as a formatting tool only.
Now that an out.2 file has been created, although filled with irrelevant data, the
user defined mode shapes need to be developed and replace the data currently in the out. 2
file. To keep the problem as general as possible, a series of EXCEL workbooks are set
up to contain each of the calculated mode shapes. Even though the mode shapes are
known, the magnitude of the modal displacements is still arbitrary. Mode 1 is simply a
rigid body plunge motion. For this mode, the spreadsheet contained a large matrix of
data containing information on the nodal displacement due to a generalized displacement
of 1 inch. Table 31 shows how the data was arranged in the spreadsheet. A similar
format was used with the other two modes.
Table 31: Spreadsheet Layout for Manual Input of a Structural Mode
Original Location New Location Nodal Displacement
Node X I y I z X' I Y' I z' ~X I ~y I IsZ
Each row in the above table contained data for every node in the solids file. The simplest
case was rigid body plunge.
coordinate.
In this case, one inch was added to each original z
39
Slightly more complicated was the determination of the rigid body pitch motion.
The magnitude of the rotation, so ~ong as the rotation was a pure rotation about the
models midchord (8 inches), was arbitrary. The modal displacement vector, however, is
sensitive to this magnitude. It makes more sense to explain this point with an example.
Consider the choice between a rigid body rotation (~2) of 1° or 10°. Two figures below
illustrate how STARS interprets the mode shapes it creates, or the user defines.
7
x
Figure 34: Rigid Body Pitch Mode Definition Example
Shown in Figure 34 is the original structural position (dashed) and the TOtated
position. Remembering from Table 31 that nodal displa.cements were specified such that
only the original and final position are known. STARS therefore interprets the entire
structural rotation as the straightline displacement from the initial position to the final
position. Figure 34 shows only the positive displacement. For an oscillatory motion,
however, both positive and negative displacements would be encountered Figure 35
shows what happens for a displacement opposite that of the defined rotational mode
shape.
40
_.

...
l·········· ..... ......
.....
.......
...•.. ~ '. ··..r
Figure 35: Structural Deformation Opposite Original Definition
As the structure rotates from its original position (dashed) to its specified deflection
(dotted & grayed), each node follows a particular vector defined by the final
displacement. As the structure rotates from this position, back through the original
position, and to the position shown in Figure 35, it follows the path defined as shown by
the vectors. One can immediately see that this structure must compress and stretch as it
cycles through its motion.
We can now see the effect of the magnitude of the specified rotational mode
shape. It makes sense then that if the rotation amount is small that any compression and
stretching can be kept to a minimum. Keep in mind also that actual displacements may
be larger than those originally specified, resulting in further contraction and expansion. It
is apparent that a compromise is needed. Structural distortion during rotation was
minimized by specifying small, 1°, rotational mode shapes, and neglecting any slight
changes in the longitudinal direction, i.e. as the object rotates, only vertical motion is
realized, translational motion is neglected. This effect is illustrated below in Figure 36.
41
....
..
....... ,
 ; ~
...... jI
Figure 36: Implemented Rotational Mode Shape for STARS
The above figure is shown at an exaggerated displacement to highlight the method used.
In actuality, the 10 rotation produces translational changes on the order of 0.001 inches.
Additionally, neglecting this small translation allows for a more general rotation angle.
For small angles, those around 80 or so, translation due to rotation can still be considered
insignificant.
Finally, specification of the control mode followed much the same procedure as
did the rotational rigid body mode. The difference being the fact that the control surface
rotated about the % chord pojnt (12 inches) as opposed to the midchord. Again, modal
displacement vectors were specified at each node, but only the nodes on the flap had nonzero
values. The same stretchi.ng/compression problems were encountered with the flap.
It was critical that the flap be modeled as accurately as possible so that any slight
deflection would be correctly interpolated to the CFD grid, hence the dense mesh in
Figure 33. The translational effects due to flap deflections were more significant than
those due to the entire wing pitching because of the relative sizes of the flap and wing.
As was done with the rigid body rotation, a 10 generalized displacement was used to
specify the motion of the flap. The EXCEL workbooks showing the nodal displacement
data are given in Appendix A2 through A4.
42

With all of the little details discussed in the previous paragraphs, it is easy to
loose sight of what has actually been taking place. Up to this point a solids mesh has
been generated, STARS has performed an undamped, freevibration analysis on this
mesh and has generated an out.2 file (for formatting purposes only), The next step is to
replace the data in the file with the natural frequencies and mode shapes that were
developed a few paragraphs back. The data in the out.2 file must be arranged in an exact
format due to a formatted read inside STARS. Inside this file, displacement and rotation
data are given for each structural node number. Displacements for each node are broken
up into x, y, and z translations and x, y, and z rotations. As opposed to entering all the
data manually, a quick FORTRAN program (Appendix A5) was written that read in the
data from Appendices A2 and A4, sorted it into the proper form and output the data into
an external file. Data from this new file is in tum manually pasted into the proper
location inside the out. 2 file. There is quite a bit of manual overhead when one chooses
to define frequencies and mode shapes that is not involved when STARS computes them.
However, time savings are realized during the latter parts of the solution when si mple
changes in mass, damping, stiffness, CG locations, etc. require the modification of a
single parameter and not the redefinition of the basic structural mode.
Throughout thi s section, there has been a reference to the CFD mesh. Before
discussing further structural requirements for the full ASE simulation, the development of
the CFD mesh must be considered.
3.2.2 BACT CFD Model Development
The development of the CFD mesh consists of several key elements. First, the
model geometry must be constructed and entered such a way that STARS can read it.
43
[ .
Next, the CFD mesh must be set such that the grid is sufficiently dense, or coarse, in the
and the solution parameters are set, the stage is set for a steadystate CFD solution.
appropriate areas. Finally, once the CFD boundary conditions are completely specified
3.2.2.1 BACT Geometry Specification in STARS
The first step in defining the CFD mesh in STARS is the specification of the lines
and surfaces that will make up both the model geometry and computational domain.
Described in the next few paragraphs is the development of two CFD meshes. The first
mesh is the is used to investigate the application of the transpiration method for a variety
of cases including steady and unsteady aeroelastic cases, steady control surface
deflections, and finaHy control surface deformations as a means of flutter suppression.
The second CFD mesh is used to compare an actual control surface deflection to one that
has been modeled using the transpiration boundary condition.
Shown below in Figure 37 are the important labels defining the lines and
surfaces of the entire computational domain.
44
~
I .
I

r·;:..]
: . f~ . ! ~"'! ~
........ !
".;
Figure 37: CFO Computational Volume Specification
In the above figure, the circles indicate the definition and specified direction of a line.
Parallelograms indicate the existence of a surface. In both cases, dasbed lines represent
lines or surfaces that would be hidden in order to facilitate the visualization of the 3D
geometry. Next, a similar procedure is employed in Figure 38 for defining lines and
surfaces on the wing itself
45
Ir.
I.
,,
.
, .
i
1
············     .. .    ~._........ .  .  .   . ..    ..  ..   
FElISA·XPLT g~ ~
y
~]
Figure 38: Wing Geometry Specification
To specify the chordwise points that define the NACA 0012 airfoil cross section, such as
lines 13 and 14, 161 cosine spaced points outline the curve of the airfoil. Cosine spacing
simply allows finer specification along the leading and trailing edges with reduced
spacing over the surface of the airfoil where there is the least curvature. Admittedly,
there is an excessively large number of points defining these curves, but any effort to
minimize any sort of modeling error was utilized. Surface 15 corresponds to the wing
tip. A rounded surface for the wing tip was included to match that of the experimental
BACT model. This rounded tip is simply a surface of revolution which is defined by 1;2
of the airfoi I section.
46
I~··_······ .._  _ ... . .__. ..   ... . .. .. .. . . "
FELlSA·XPLT ld~m
o
<·.~·.··1~,,·:"'1~
". . .
.............. ;.
Figure 39: CloseUp of Deflected Flap Geometry Definition
Figure 39 demonstrates the additional lines and surfaces needed to define the
wing geometry for a deflected control surface. For clarity, only lines and surfaces that
were modified or added were included in this fi.gure.
Addition of the control surface causes several difficulties due to the changing
intersection points between the control surface and the wing. With any slight change in
the deflection angle, intersection points must be recalculated and the STARS data file
modified. This manual remeshing concept is, as one would expect, time consuming
The time it takes to go from one deflection angle to another is on the order of 2 to 3
hours. That is just the time it takes to modifY the wing geometry. Changes also must be
made to the file that contains information about grid density and element source location
47
since surfaces are being displaced from their original positions. Additional time must
also be spent regenerating the tetrahedral mesh throughout the entire computational
domain. This process itself, can take a couple more hours. All in all, the time it takes to
go from one deflection angle to another can take on the order of 5 or 6 hours to
completely redefine the geometry and regenerate the computational volume. The
complete data file is given in Appendix B] .
The procedure described above serves as a very good basis for the use of the
transpiration method to model these types of discontinuous deflections. In an
environment where it is desired to obtain results for a number of control surface
deflections, one could easily make the simple modification to the scaling factor that
describes the control surface deflection angle. For example, a 0° deflection is equivalent
to saying "Zero times the generalized displacement of 1°." Similarly, for a 10°
deflection, it equivalent to saying, "Ten times the generalized displacement of 10
" It is
important to note here that a positive control surface deflection angle corresponds to a
downward deflection ofthe flap.
As in the SOLIDS definition, a series of EXCEL workbooks was set up in order
to facilitate the assembly of the data file STARS uses to create the surface front. The
spread sheet is set up in such a way as to automatically redefine each surface and line
definition for any symmetric 4digit NACA series airfoil crosssection. Due to the
number of reference points defining the airfoil crosssection, the data file is nearly 6000
lines long. One can immediately appreciate the use of the automatic file generator for
such a large number of points. For the case of the actual control surface deflection,
however, the data file generation cannot be done automatically. With each deflection
48
angle, the intersectton points discussed earlier change, so one must go through and
compute the intersection points and redefine lines 24, 26, 33, 34, 36, 37, 39 and 40.
3.2.2.2 BACT GIid Specification in STARS
To this point we now have only the lines and surfaces that define the CFD
geometry. Next, we need to specify the location and density of the tetrahedral elements
that will define each surface and the internal volume. STARS allows one to specify
point, line, or triangular sources. These sources can be thought of as sources of
tetrahedral elements. Based on the specifications, tetrahedral elements will originate
from the point, line, or triangle at a specified density and taper off toward larger elements
based on another specification. For the BACT wing, line sources were placed along the
leading and trailing edges of the wing, the upper and lower surface locations that
correspond to the beginning of the control surface, and along the wing tip. An
arrangement of triangular sources lie under the surfaces of the wing and control surface.
Arriving at an optimal grid density is an iterative process. One simply begins
with a grid that seems right and iterates based on the mesh observed. With this file
specified, STARS is able to assemble the mesh for each surface which can then be
viewed to get a visual sense of the grid density. The resulting mesh for the BACT wing
can be seen in Figure 310. This figure contains four different views of the surface mesh.
The mesh is dense where one would expect high flow gradients, and less dense where
there the flow gradients are not as sharp. Of course, it makes sense to have a more dense
grid at the leading and trailing edges, at the wing tip and over the control surface region
but the grid density over the upper surface of the wing seems overly dense at first glance.
This is explained by the simple fact that the BACT wing was tested at transonic Mach
49
numbers. At a Mach number greater than Mach 0.77, transonic shocks begin to appear
on the surface of the wing. As the flap rotates up and down, these shocks also translate
across the surface of the wing. During flutter, as the wing pitches and plunges, the
location of the shock changes once again. In a full aeroservoelastic simulation where the
wing experiences each of the cases mentioned above, the location of the transonic shock
can manifest itself at almost any chordwise location. Also, for this relatively low aspect
ratio wing, one would expect the threedimensional effects to be significant.
Figure 310: Views Showing Tetrahedral Surface Mesh on the BACT Wing
Therefore, in order to accurately capture the full threedimensionality of the flow and the
location of the transonic shocks, the grid density over the entire surface of the wi ng must
so
be kept sufficiently dense. The file containing the specifications on the location and
dens~ty of the tetrahedral sources is given in Appendix B2.
What has been constructed thus far are the wing and flow domain lines and
surfaces and the surface discretization for each surface. What is lacking now are the
three dimensional tetrahedra that will constitute the rest of the computational domain.
What we were able to see in Figure 310 was the grid density on the wing and wall.
Common sense dictates that the more tetrahedra one has in the flow domain, the longer
the solution will take to converge. There is, therefore, a tradeoff between a sufficiently
dense grid and solution time. The authors of the mesh generation code recommend
Equation (33) as an approximation of the number of mesh points as a function of the
number of surface points [peir6, Peraire, and Morgan, 1993].
Where N~ = Number of Mesh Nodes
C = Empirical Constant (1.62)
N; = Number of Surface Nodes
n = Empirical Constant (1.15)
(33)
Table 32 shows a the number of surface nodes and a comparison between the number of
mesh nodes resulting from running the volume generator for the BACT model, and the
suggested value from (33). Additionally, the number of tetrahedra in the computational
domain is 342,469.
Table 32: Actual and Suggested Number of Mesh Nodes for STARS Volume
51

The 63,902 mesh nodes compares reasonably weB with the 55,778 nodes predicted by
Equation (33).
3.22.3 BACT Boundary Condition Specification in STARS
From 3.1.2 we see that the next step is to run the SETBND routine to define the
boundary conditions for lines and surfaces. This routine uses the file, found in Appendix
B3, to specifify walls, farfields, symmetry planes, singularity lines, etc. STARS uses
this data to assign the proper CFD boundary conditions on the nodes adjacent to the
specified elements. For the BACT wing, the back wall and all of the wing surfaces are
defined as walls. The remaining surfaces are defined with/arfield boundary conditions.
Lines along the trailing edge are defined as singularity elements. A singularity line
simply defines a region in the CFD model which does not have a well defined normal,
such as the trailing edge of th.e wing, where the upper and lower surfaces end at a sharp
point, there is no way to specify a single normal. Ignoring singularities can result in
abnormally high flow gradients that tend to washout the true flow physics.
The last thing that needs to be done is to specify constants that the flow solver
will use throughout the solution. This is done using two files. The first is the CONU file.
This file specifies the number of timesteps to run, the number of innerloops to run at
each time step and a host of other parameters. This file is given In Appendix B4 so only
those parameters that are of key interest to running a steady solution for the BACT case
are discussed.
52

3.2.2.4 Effect of the Dissipation Parameters in STARS
Making use of the inviscid flow assumption can be problematic in the transonic
flow regime. Here, transonic shocks on the surface of a wing tend to be weak. With an
Euler solver, these shocks tend to be predicted later and more sharply than shown with
experimental data. STARS allows the variation of a few control parameters that
introduce dissipation into the numerical solution. Changing the values of diss(l) and
diss(2) in the file BACT.CONU, given in Appendix B4, had a very significant impact on
the pressure distribution prediction. From their default value of 1, the constants were
eventually modified to their current value of 3.5. Figure 311 shows the predicted
pressure contours, with and without modified dissipation constants, compared to those
I .•"
obtained through experiment.
1.0
0.6
o.S . ........·c)....I. '' \
~ ~.
0.4 •• .j11
0.' ~ , • .,..".....:*., G0.0 a•
• 0.2' , .... ••
0.4
0.6 .
0.8
• Experimental (M=O.82)
ASTARS  DiSSIpation
• STARS  No DissipatIOn
0.5 0.6 0.7 0.8 0.9 1.0
xlc
01 0.2 0.3 0.4
1.0 +~tI.~.~+~I~~+~~t~+"~~+~""'r''~+~""i
0.0
Figure 311: Effect of Dissipation Constants on Cp in STARS CFD Solution
53
As the figure illustrates, the predicted transonic shock without dissipation is predicted aft
of the actual shock and is more sharp in nature. Including dissipation allows for very
good agreement between experiment the STARS prediction.
Determination of the best value of the dissipation constants was an iterative
process. For the range ofMach numbers at which the BACT wing was tested, the highest
value of dissipation that did not cause the solution to go unstable was ~3.5. Dissipation
was not noted to improve the solution convergence time, which is discussed in more
detail next.
3.2.2.5 SteadyState Solution Convergence Criteria
As the steady solution starts out from a given freestream Mach number, the
resulting flowfield about the geometry evolves through time. As the solution progresses
STARS outputs residual values. These residuals are an indication of how much a flow
parameter, such as density and velocity have changed since the beginning of the solution.
Typically, once the residuals "become small enough", the solution is said to have
converged. What was discovered with the BACT wing, however, is that the residuals
were not necessarily the best indicators of convergence. The item that ended up being the
most convenient indicator of solution convergence was the maximum Mach number. The
judgment of when the residuals were low enough was too subjective. The maximum
Mach number gives a more objective view of solution convergence. To further make thi.s
point, a comparison between the two methods is given in Figure 312.
54
,0<. ,oo, ......ILLL._~__~~
~ 011. 1OIll I UIIlI
Jm • ....
11.ID.l ~=
·....,
Figure 312: Solution Convergence Using Residuals and Maximum Mach Number
The picture on the left in the above figure shows how slowly the residual drops for a
given case. Even on a log scale, there is no definite solution convergence. The picture
on the right, however, clearly shows that the maximum mach number converges to one
particular value.
3.2.3 BACT Uncertainty Estimation
Before the development of the aeroelastic and aeroservoelastic models are
developed, one must consider the experimental uncertainty present in the BACT model.
As with any experimental measurement, we expect to see a certain amount of
experimental uncertainty_ These experimental uncertainties, unfortunately, were not
quantified for the BACT model. In an effort to determine estimates for these
unceriainties, communication with Mr. Robert C. Scott and Mr. Martin R. Waszak of the
NASA Langley Research Center, provided valuable insight into the uncertainty of the
measurement techniques.
Since the BACT wing is considered rigid, aU of the stiffness terms arrive from the
use of the pitchandplungeapparatus (PAPA) [Farmer, 1982] _ The wing is reportedly
mounted on the PAPA such that the elastic axis is coincident with the geometric center of
...
55
the PAPA mount. In the next few paragraphs, estimates in uncertainty are given for the
determination of structural mass, stiffness, and damping characteristics, the location of
the center of gravity relative to the elastic axis, and determination of dynamic pressure at
flutter.
First, estimates In the uncertainty involved in the determination of structural
stiffness and damping is covered. For a twodegreeoffreedom model, the stiffness
terms of primary concern are the plunge and pitch stiffness. To measure the plunge
stiffness, weights were attached at a location corresponding to the wing's midchord.
Stiffness was then determined simply by dividing the additional weight by the resulting
deflection. Similarly for the pitch stiffness, a known torque was applied about the wing's
midchord. This known torque was divided by the resulting angular displacement in
order to determine the pitch stiffness. Structural damping was determined by exciting the
structure in either pitch or plunge and measuring the decay in the freeresponse
Generalized mass of the pitch and plunge modes was determined from the
resonant invacuo natural frequencies. The resonant frequencies were determined by
exciting the structure in either pitch or plunge and measuring the number of cycles in a
fixed time. Knowing that the natural frequency, stiffness, and mass are related by, (34),
one can calculate the generalized mass from the measured stiffness and natural
frequency. The resulting measurements from the above tests are summarized in Table 3
3. Literature only reports those values in test # 3, the author appreciates the additional
data from Mr. Waszak.
(34)
56
Table 33: Experimental Measurements in Structural Parameters
2 2637 2964 0.00150.0018 3.360 5.302 6.03 2.70
3 2686 3000 0.0014 0.0010 3.344 5.208 6.08 2.80
Recall that the Benchmark Models Program at NASA Langley involved tests on
both a NACA 0012 wing as well as the BACT wing, both tested and mounted on the
PAPA with the wing's midchord nearly coincident with the elastic axis. The two wings
had the same chord, span, airfoil crosssections, and experimental instrumentation layout.
The only external differences that exist are small geometric defects, and the presence of
three control surfaces. Internally, a portion of the material had to be removed for the
installation of actuators etc. Despite the material removed to add the actuators and
spoilers and separate the trailing edge control surface, structural characteristics are very
similar between the two. Rivera and others report the values shown in Table 34 for the
structural properties of the NACA 0012 wing and PAPA mount [Rivera, et al. 1991 &
1992].
Table 34: Experimental Measurements in Structural Parameters
i'~:~~~~jl:'I' ;lJ1:·~~ijl!:·I·1'·rWll~ijl!·'!l!;:!fi!~:;lli·i·;':~·ii;"·~~f,:I·,!11:il;i!I.~~jl!i··,l·!;:l:';'~~::":;!! ;j.·il~fJ~jf1f.fli;!J:M~~&~t~~!!:.
3/92 2659 2897 0.0024 0.0024 3.36 5.20 5.966 2.714
7/91 269722854.6 0.0034 0.0016 3.40 5.18 5.910 2.7695
Recall that values shown in Table 33 were obtained from the BACT wing and PAPA
mount. Comparing these values, we see that the tables are very similar. This would
seem to indicate that the physical differences between the model should be essentially
negligible.
57
....
Next, experimental uncertainty in the determination of the center of gravity (CG)
relative to the elastic axis (EA) proves to have a very significant effect on the prediction
of the flutter boundary. In the literature, the CG's location relative to the elastic axis was
reported, at best, to be nearly coincident with the midchord of the wing. Waszak reports
the value of the inertial coupling between the pitch and plunge modes (Sh.a) as being
0.0142 slug·ft. Using (35) below, we can estimate the relative location of the CG to the
EA.
S =m·x Cl cg (35)
Using the value of Sa. reported by Waszak (1996), and the mass, we calculate that the
distance from the EA to the CG is 0.028 inches. After communication with MT. Waszak,
he stated that his reported value of 0.0142, which was well within experimental
uncertainty, had to be used to account for the slight difference between his computational
model and experimental data.
The presumption that the CG and EA are coincident comes from qualitative
observations made during testing at the NASA Langley Research Center When
measuring the plunge stiffness, weight was applied at the midchord of the wing. During
these tests, there was no reported difference in the displacements of the leading and
trailing edges indicating the absence of static coupling. Similarly, during measurement of
the pitch stiffness, no static coupling was observed. In order to excite the pitch and
plunge frequencies, the BACT wing and PAPA mount were excited by an initial
deformation that was suddenly removed to allow the structure to vibrate freely. This was
done in a manner similar to the static loading, at the midchord for plunge and about the
midchord for pitch. The free vibration of the model excited in this way showed very
58
little pitch motion when excited in plunge in very little plunge motion when excited in
pitch. This, of course, implies that the CG must be very close, if not coincident with the
midchord. As an approximation, the relative location of the CG and the EA was
estimated to be no more than 0.1 inches, or 0.625% of the chord.
From the above data, it is possible to construct a simple model from which we
could quickly evaluate the importance of some of the above parameters. Using a
simplified aerodynamic model and solving the equations of motion using the pmethod,
we can quickly solve for the divergence speed. Since the flutter speed can be solved for
explicitly in the pmethod, parametric studies can be done very quickly. Shown below in
Figure 313 we see the effect of three different parameters on the flutter prediction.
100% 0.75% 0.50'/. 0.25% 0.00% 0.25% 0.50% 0.75% 1.00% 1.25% 1.50%
+ Plunge Stiffness
...... Pitch Stiffness
Ar CG Location Relative to Elastic Axis
'I· I' ., I , I
50%
2.5%
00'/0
2.5%
t;
50% {
'" ,
" ·7.5%
~
.S !O.(l%
~Jl
§ 12.5% :
= U 150% ..
':f? Co
17.5%
20.0%
225%
250% I
1.50% 1.25%
% Change in Parametl~r
Figure 3 13: Effect ofKh, Ka, and xcg Location on Flutter Prediciton
In the above figure, the x axis represents small deviations from nominal values for plunge
and pitch stiffness, Kh and Ka , and xcg , which is a measure of the distance from the elastic
59

axis to the center of gravity, measured positive aft. As is shown, small changes in both
plunge and pitch stiffness effect little change in the flutter prediction, ±2.5%. Small
changes in xcg , however, influence the flutter prediction significantly. Using the above
model, changes in xcg on the order of 1% can change the flutter prediction by over 20%.
Now, having the BACT's CG location specified as nearly coincident with the
elastic axis introduces a slight difficulty in flutter prediction using STARS. For
comparison with experimental data, small variations in each of these parameters can add
up to large differences in flutter prediction. In addition to all of these differences, there is
still the matter of determining the actual flutter point. Looking at time traces of
experimental data, it is often difficult to tell exactly when the system is going unstable.
Mr. Waszak estimated that the dynamic pressures that defined the flutter boundary were
measured to within ±2 Ib/ft2
.
Knowledge of these and other uncertainties is fundamental to appreciating the
degree to which the computational model can approximate the experimental data. When
developing an aeroelastic model, we must assume that the wing is exactly rectangular,
perfectly symmetric, its crosssection exactly matches that of a NACA 0012 airfoi I, its
mass, damping, and stiffness, and the coupling between each, is known precisely, etc.
One can quickly appreciate the amount of tolerance buildup that is present in the
experimental data. These small, relatively lwknown, differences translate into a lot of
finetuning of the computational model. In work presented by Waszak, flutter prediction
within 7% of experimental data was considered" .. pretty good ... " [Waszak, 1998).
60

3.2.4 BACT Aeroelastic/Aeroservoelastic Development
From 3.2.1, 3.2.2, and 3.2.3, we now have a solids model with three specified
mode shapes, a CFD model, and an appreciation of the experimental uncertainty involved
in the aeroelastic data. As developed previously, aeroelasticity is the coupled response of
the two aforementioned models. From section 3.2.2 we have the capability of producing
a steady CFD solution from which to begin an unsteady simulation. Next, the mode
shapes specified in the SOLIDS module for the finite element structural mesh must be
interpolated to the CFD mesh. Using the orthogonal property of the natural mode shapes,
a superposition of these natural mode shapes can be used to represent an arbitrary
structural deformation.
As mentioned previously, if the natural frequencies, mode shapes and other
structural properties are known beforehand, they may be entered manually into STARS.
Before the interpolation begins, STARS must know which surfaces represent moving
boundaries. This is done with information contained inside the ftle BACT.SCALARS,
given in Appendix B5. This file contains a variety of other parameters of interest to the
unsteady solution, but those are not of particular interest and will not be covered. When
the interpolation from the SOLIDS mesh to the CFD mesh occurs it creates an AJUU YS
file. This file contains information regarding the natural frequencies for each mode
shape, the generalized mass, stiffness, and damping matrices and nodal displacements for
each CFO node which represent nodal displacements on surfaces in the CFO mesh for
each mode shape. The BACT case has 3 modes. plunge, pitch, and a control mode that
represents the moving control surface. The mass, stiffness, and damping matrices are
therefore 3x3 matrices. There are three sets of nodal displacement, or AERO vectors,
61
one for each mode that specify the generalized disptacement of each CFD node for each
mode shape. Only the top portion of the BACT.ARRAYS file is given in Appendix B6
because the file is over 26,000 lines long.
In STARS, the mass, damping and stiffness matrices were manually entered into
the BACT.ARRAYS file such that they matched those reported in test #3 for the BACT
wing. Since geometric data was entered into stars in units of inches, units of mass are in
slinches as opposed to slugs. Where a slug has dimensions of Ibfs2/ft, a slinch has
dimensions of Ibf·s2Iin. The conversion is, therefore, I slinch = 12 slug. Observing the
plunge equation we encounter no dimensional conflict within STARS. Noting the
moment equation, (36)(313), we see the possibility for a slight discrepancy. Beginning
with the general equation for the moment in (36) we see the following.
Id+f3a+Ka=M (36)
The moment is simply the integral of The moment is simply the integral of the pressure
times the mode shape, so substituting this into (36) we arrive at (3 7).
fii + f3a + Ka =fp¢x.ix
In STARS, however, we have the following definition, shown in (38).
fp¢dx =AI/ao =a ofp¢dx
(3.7)
(3.8)
(3.9)
Where ao, is the generalized pitch displacement of 1°. Rearranging the equations, we get
STARS definition of the pitch moment in (39).
aolii + a of3a + aoKa =fp¢xix
62

Solution of (39) assumes a in units of radians, so using (310) we must convert the
angular displacements and velocities displacements into dimensional form consistent
with the generalized displacement.
We can now substitute this relation back into (39) and obtain (311).
aJ(~aJ + aofJ(~aJ + aOK(~aJ =fp¢x:ix 180 180 180
(3.10)
(311 )
Now, units are consistent on both the left and righthand sides and can be arranged into a
more convenient form, shown in (312).
(ao~IJa + (ao~f3Ja +(ao ~KJa = fp¢x:ix 180 180 180
(312)
Remembering that the general ized pitch displacement was I° or nl180 radians we can go
ahead and multiply the generalized displacement by the n/l80 factor and obtain (313).
(~l)a +(~fJ)a +(~K)a =fpr/xlx 180 180 180
(313 )
Since generalized displacements for the wing and flap are specified as 1°, parameter entry
into the system matrices within STARS requires a premultiplication by n2!l802 Note
that this problem was not encountered for the plunge degree of freedom since both the
generalized displacement and modeshape are in inches.
Shown below in (314) is the mass matrix that is entered into the BACTARRAYS
file. Notice that rows 2 and 3 are premultiplied by the n2IJ 802 scaling factor.
63

Where:
m Sh,a ShS
2 2 2
~S ~I ~S
1802 a 1802 a 1802 a.S
? ~ n 2
~S n
S 1 1802 h,S 1802 a.S 1802 /j
m  Generalized Mass (Plunge)
Ja  Generalized Mass (pitch)
ls Generalized Mass (Control Surface)
Sh,a  PlungePitch Inertial Coupling Term
Sh,o PlungeControl Surface Inertial Coupling Term
Sa,o PitchControl Surface Inertial Coupling Term
(314)
Similarly, the damping matrix is shown in (315):
gh 0 0
0
n 2
0
180 2 fL
?
0 0
n
] 802 g J
(315)
Where gh  Generalized Plunge Damping
ga  Generalized Pitch Damping
go  Generalized Control Surface Damping
In the above relationships, g is defined to be M2 t; 0Jn, where M, S, and 0Jn are the
appropriate generalized mass, damping, and natural frequency.
....
o
o
o
o
64
o
o (316)

As mentioned earlier, the uncertainty present in the BACT afIects the values that are
entered into (314) to (316), From the sensitivity study, \\Fe saw that the most sensitive
uncertainty exists in the specification of the pitchplunge coupling term So. since this
related directly to xcg as discussed in the previous section, Final system matrices were
obtained after finetuning the parameters and are given in Table 35.
enera 17.e ass atn.x
0.50667 0.0506667 O,00n8
o154339x 104 0.0102350 o57390x 10')
0.877298x 10'0 0.57390x 10'5 0.40134x 10. 1
Table 35: BACT Model Parameters in STARS
G f dM M
GenerarlZetIDammnf! M.alnx
0029819 0.0 0.0
0,0 0,66985 x] 0'3 00
00 0.0 0,0
Generalized Stiffness Matrix
0.223833 xl 03 0.0 0.0
00 o I09500x 10 00
00 0.0 o 1096623x 104
3,2.4.1 TimeStep Definition in STARS
The timestep used in STARS is computed using (317) where the parametersfreq
and flslpe are defined in the CONU file, M is the freestream Mach number and a is the
sonic velocity
d/ = 21f
freq "nstpe . (A1 .a)
(317)
Until recently, there has been no prescribed method of determining the timestep, A
generally accepted ruleofthumb was to make certain that at a single period of oscillation
at the highest frequency was made LIp of at least 3040 timesteps For supersonic cases,
65
>
this seems to work just fine. In subsonic flow, however, wake effects are propagated
throughout the entire computational domain. More recently, fi'eq, is defined to be simi lar
to the highest natural frequency in the SOLIDS model. The parameter nstpe can then be
thought of as the number of time steps per period of oscillation. One can also look at this
another way. In STARS, the default value of nstpe is I. Instead of letting/req represent
the highest frequency, it may be arbitrarily set such that one obtains an equivalent timestep
within STARS. Either method works equally well, but lettingfreq represent a true
frequency and increasing the number of steps per period (nstpe), makes more intuitive
sense.
for proper flow dynamics, the user is concerned with the number of inner CFD
iterations per time step (ncycl), and the length of the timestep. A recent investigation in
STARS with a simple NACA 0012 airfoil provided valuable insight into the relative
importance of ncycl and dt. The study was done using the problem of a suddenly
accelerated wing in subsonic flow (Wagner Problem) wh~re the lifL and drag arc timeevolving
parameters. While changing the parameters t1cycl and dr, plots of the changing
lift were obtained and plotted vs. a nondimensional time parameter. Each of the
following plots were obtained for an impulsively stalied ACA 0012 airfoil at Mach 0.3
at a.=5°
As Figure 314 demonstrates, for a given value of ncycl the time varying lift is
highly dependent on the size of the timestep The disadvantage of going with a small
time step is, of course, the fact that it will require additional computational time to run an
equivalent job which incorporates the larger time step.
66
zC"f

Another option is to keep the timestep the same, but let the CFD solver perform
more iterations at each timestep. This case is demonstrated in Figure 315 where we see
the effect of time step on differing values of ncyd Shown here is, again, a high degree
of sensitivity to the size of the timestep for a given value of ncycl. For the timestep of
0.1 in the upper plot, changing the value of ncycl shows a definite etfect. The lower plot
shows that for a much smaller time step, 0.025, the plots of CI vs. t* are virtually
identical despite the fact that values of ncycl differ by a factor of 4. The conclusion,
therefore, is that the importance of a small timestep outweighs the imp0l1ance of
increasing the number of CFD iterations per timestep.
67

0.8 ..
u \ 0.6
\.
0.4 . .,.,•
0.2·
"Exacl"
&NCYCl.=IO, <1,"=0.1
NCYCl.=IO, 1Il"=O.02S
AsymplO1C

,.\ I b ,.
u
14
L2
0.8
06
0.2
..
~'Exacl"
<JNCYC!=40, <11"=0.2
+rvcrCl.=40, <11"=0.1
NCYCL=40,61"9).025
/\SymptI11c
0' .~~ ~~,''''"'  ._. , ,"'
02 00 08 ,. 12 JA 16 \8
.....
Figure 314 Effect ofTirneStep at Two Different Values of ncycl on the Lift Evolution
for an Impulsively Started NACA 00] 2 Airfoil at Mach 0.3, a=5°
68
J.'

J'
~..
[)
~.6
0.4
0.2 .
"0 0.2 0.4 O' 0.'
,.J
, 2 .
"Exact"
NcrCL~IO,JiI''''ll1
+NCrCLdlO, LlI·~O.1
I>NCrCL~40,Jit""'ll I
.4..symplOle
t'
......... NCYCL=40, M""'l1015
 "cxaC{"
~ NCYCL=IO, ~I·~O 015
08
c ~
"' \ "~'",""
04L + I ~Oo'O I ............................
02
o ~ ~~.,_. ,,~, ._. , . ' .. ~_._. ,
02 0·' 0.' o.
,*
1.1 ,,' 16 I,j
Figure 315' Effect of ncycl with at Two Different TimeSteps on the Lift Evolution for
an Impulsively Started NACA 0012 Airfoil at Mach 0.3, a=5°
To begin the ASE simulation, we must have first generated the ARRA YS file and
completed a steady state CFD solution at the reference Mach number Once the
parameters are set in the SCALARS file, and the CONU file is configured properly, an
ASE solution may be started The length of the solution is determined primarily by
parameters in the CONU flle. There are a lot of parameters set in this file, but for the
69

ASE simulation, we are primarily concerned with the values ofnstep and ncycl. The total
number of timesteps is specified with nstep The number of inner CFD iterations per
time step is specified with ncycl. For instance, with nstep = sao, and ncyd = 40, the ASE
simulation would last for a total of 500 outer time steps. At each inner timestep, 40
CFD iterations are allowed for the computation of the new aerodynamic forces. All
together, these parameters specify that 500x40 CFO iterations.
3.2.4.2 Modal Identification Technique
With each CFD iteration taking on the order of 30 seconds, we quickly see how
timeconsuming these ASE simulations are. For the BACT, the Ilstep and f/LYcI were
generally 5000 and 40. Assuming 30 seconds per CFO step and doing the math, we
estimate that an EULER solution for a single transient may take on the order of 69 days
on an IBM RS6000 38T The general procedure required that the solution be monitored
and when the timehistories looked to be going unstable, assume that is the flutter point
and kill the solution. The nature of the BACT system makes this method impractical
Observe the following figure which is a portion of an actual timehistory obtained hom
STARS.
70

0.6
/
/
I._~.  1. ~j. I !
0.4 0.5
 Mode I (plunge)
__Mode 2 (pitch)
r //'\\/ //''
I
I ' ....~ ~,/ /\ ,/'
, \. i , J '\""'/
 1.5 I:''r,,~I..LJ .. ~',~LL;' , _, _ ,
0.0 0.\ 0.2 0.3
Time
1.5
1.0
~ 0.5
"'0 .a :.:::l 0.0
~ 0.5
1.0
Figure 316: Abbreviated TimeHistory of BACT Wing in STARS
From Figure 316, it appears that the solution is going unstable. Typically, that would
have been considered good enough but allow the solution to continue Jor the full 5000
timesteps and we obtain the following timehistory.
Mode I (Plunge)
__Mode 2 (Pitch)
1.5 2.0
Time
2.5 3.0 3.5
r'
I ( ,J
, I! \; \,
! I
I
, i
l\ j
1.0
)!
II,,
0.5
I ,
i r I'i
1.5
1.0
Il.l 0.5
"'0
.E :.:::l 0.0
Q.,
~ 0.5
1.0
1.5
0.0
Figure 317: Complete SOOOStep TimeHistory of BACT Wing in STARS
Allowing the solution to continue for the full 5000 steps, Figure 317 shows that mode 1
exhibits a slight amount of dampedbeating whiJe mode 2 is lightly damped, therefore not
...
7I
yet at the flutter point. Beyond visual interpretation, a modal identification technique
provides damping characteristics for each structural mode [Eckhart, 1998]. Given a timehistory
from STARS, this tool provides the user with both a damping frequency and
damping factor.
3.2.4.3 System Identification Technique
.....
For each flutter point, one must generally take a trialanderror approach in the
determination of the flutter point. Trials are made until dynamic pressures on each side
of the flutter point are obtained. This is, of course, very time consuming. The
determination of a complete flutter boundary for a problem of this type could easily take
several months Recent work by Cowan, allows the use of a system identification
procedure to model the coupled structural/CFD system (Cowan, 1998]. This has the
signifi.cant benefit of accelerating the time required for a full ASE simulation. Essentially
eliminating the CFD solver, which makes up the vast majority of time during a coupled
simulation, and replacing it with an algebraic transfer function reduces ASE runtimes
from days and months to minutes. For the same 5000 step solution described previously,
an ASE simulation is obtained in about 5 minutes, as opposed to 69 days.
This system identification procedure is currently implemented into the STARS
and provides a very accurate prediction of the full Euler solution To model the system,
each mode is displaced from an initially steadystate CFD solution through a known input
referred to as a multistep Forcing the CFD model with these known modal inputs
during an Euler solution allows STARS to compute the aerodynamic forces due to these
known inputs. The system identification procedure then constructs an ARMA model
based on the known inputs and resulting outputs. Once the system is modeled, the Euler
72
aerodynamic solver IS essentially "replaced" with a much faster system of algebraic
equations.
The multistep sequence used on the BACT wing is given in Figure 318 and
specified with parameters in the SCALARS file. The duration of the multistep is
determined by the following equation: 5 + isize(4nr +3), where isize is the magnitude of
the multistep and nr are the number of modes to be excited. For the BACT, isize and nr
were generally set as 10 and 3, respectively resulting in a duration of 155 timesteps. The
actual CFD solution extended to 240 timesteps to insure that all of the aerodynamics
have enough time to come to steadystate values.
..
'J
,,' ~
.: \,
."
.__Mode I (Plunge)
__Mode 2 (Pilch)
....... Mode 3 (Control Surface)
.:.,
" ,
i
r
I
I ,,
0.005 i ;
I !
! I
0.030
0025
 0.020
Q=="l
E
Q"l u O.OISi
CIl I
C. .~ 0010 0
0.000 ,I,1_''_'''~_~
0.00 005 0.10
Time
0.15 0.20
0.15 0.20
. Mode I (Plunge)
__Mode 2 (Pilch)
....... Mode 3 (Control Surrace)
0.10
Time
I I
0.05
.1.000 ",
0.900;
0.800 i
~. II
0.700 I '
.c' 0.600 : I
.tj 0.500
.2 > 0400.
0.300
0.200
0.100
0000 '__l~
000
Figure 318: MultiStep Sequence for the BACT Wing (3Modes)
73
Notice that each mode is forced to undergo a displacement resulting from a specified
velocity. These displacements and velocities, through the transpiration method, are
implemented as unsteady boundary conditions in the CFD flow solver. The resulting
aerodynamic forces and moments resulting from this sequence of events are then
modeled. The extent to which these models actually fit the data is described in more
detail in Chapter 4, but Figure 3] 9 shows a comparison of the actual and modeled
response to the multistep shown previously.
o Model
Euler
1.4 _.
16 I' , I'" 1"" I "
0.00 0.02 0,04 0.06 0.08 0,10 0.12 0,14
0.2 .
0,0 
,.. 0.2"
~ 0.4~
.::: 0.6  d. 08
.... 1.0~
12
Time
o Model
 r':uler
0.04 0.06 O,Og (j,IO 0.12 0,14
Time
006
~ 0.04
~ 0,02
0
~ a,oo
"5 002 ~
N 0,04
~
~
0,08 
0,00 0.02
Figure 319, Modeled and Actual Response to MultiStep Input
Similar results are obtained for all other Mach numbers under consideration and are given
in Chapter 4 Further validation and references are found in the original work by Cowan
74
....
3.2.4.4 Control Law Development
The final objective of the current work is to use the trailing edge control surface
on the BACT wing as a means of flutter suppression In a paper by Waszak, the BACT
wing is modeled at Mach 0.77 in a MATLAB program [Waszak, 199697] The program
developed essentially provides the user with a statespace representation of the
BACT/PAPA system at a userdefined q at Mach 0.77. Using only a pOl1ion of the
program, models of the BACT/PAPA system were obtained at three different dynamic
pressures: a little below, close to, and beyond the flutter point. The resulting statespace
system was then condensed down into a single t,Jroup in SIMULINK. Shown in Figure 3
20 is the complete model developed with its core element being the qc1ependant statespace
model obtained from the MATLAB program.
Pitch
Plunge
Flap Angie
Output Da\
1~
Mux
Flap
Response
dis13
StateSpace Models
offlon2
111+~~
+ ~ Disturbances
System
Inputs
01+.1
Constant
q=158.62 q=12978
0}
Flap Def
~~Pi
~
on/off Product off/on !FI::~::ut
controller
Control Surface
Inputs
Figure 320 MATLAB/SIMUUNK~MocleJ of BACT with Control
...
75

While looking complicated, this is a relatively simple block diagram of the entire system.
The entire diagram basically consists of four pal1s: state inputs and outputs, disturbance
inputs, controllable and control surface inputs, and a means of viewing the output. In the
center of the diagram is the BACTfPAPA system with its associated inputs and outputs.
In the upper right are the disturbance inputs from which one may disturb pitch, plunge,
and a host of other parameters. The upper left of the figure is essentially the control
p0l1ion of the diagram, where the deflection of the control surface is controlled through
simple P, PI, PO, or pro control based on pitch and/or plunge rates or displacements.
Left of center are separate control surface inputs. If control is turned off, arbitrary
control surface displacements, sine waves etc., can be input into the system. The lower
righthandside of the figure contains blocks that display pitch, plunge, and control
surface deflection as a function of time. This too] was used to gain an understanding of
the effectiveness of different control laws before their implementation into STARS
Shown below in (3 J 8) and (319) are the equations for lift and moment of the
BACTIPAPA system employed in the model shown above.
M =qScc.1 =qSC[C'1 + C' I a +C\I 6 + z: .{C~I a +CH e+ ('\f. 8)1 (319)
.\0 I" f, 11; (1 • ,It 2(I ~  ,~ . ({ , ,~
 {,
Static aerodynamic parameters were obtained from experimental data and previous windtunnel
experiments, force and moment data at various angles of attack and control surface
positions were used to compute most of the stability and control derivatives, while
dynamic derivatives were obtained from computational analysis. Parameters unknown or
unavailable were assumed to be zero.
76

Though simplified through modeling assumptions, the model proved to be a very
useful tool in obtaining quick qualitative data regarding flutter suppress.ion using the
trailing edge control surface. Despite the quality of this data, the majority of ASE
simulation was conducted in STARS since the modeling simplifications are not a limiting
factor. Of particular interest are the additional aerodynamic mass and damping terms that
result from the plunging motion of the wing and the effect on lift and moment due to the
rate at which the control surface deflects.
In general, any control law will have as its output a desired flap position. Waszak
reports the control surface actuator's transfer function as (320) where k (1.02 deg/deg) is
the actuator gain, (56) is the damping ratio, OJ (rad/sec) is the natural frequency, Os is
the desired control surface deflection, and 0 is the actual resulti ng deflection
o kOJ
==,
OS S2 +2r;OJs+OJ 2
(320)
For our purpose, however, it is more convenient to move from the fi"equency domain
back to the time domain. The corresponding differential equation is shown in (321)
(321 )
To begin putting the above equation into statespace format, we'll make the following
substitutions: XI =0 and x2 =0 = XI' Taking derivatives of these equations results in
the following. :(1 =5 =x 2 and x2 =0 = x\ Using these relationships, we rewrite (3
21) in the following form x2 +2r;OJx2 +OJ 2
X 1 =kOJ 2o,. We now have two firstorder
differential equations which we can write in statespace format, see (322)
(322)
....
77

To actually implement these equations into STARS, the timederivatives are replaced by
the following relationships shown in (323) and (324) where n is the current value, and
n1 represents past values.
(323)
(324)
Solving each for the current values of x, and X2 yields (325) and (326)
(325)
(326)
Up till this point, the desired control surface position has been arbitrary For our purpose,
the desired flap angle will be a function of plunge displacement and velocity, and angular
displacement and velocity. The resulting control law is shown in (327), where the gains
K, are not necessarily absolute. Given the range of Mach numbers, it is assumed that
some sort of gainscheduling, based on both Mach number and dynamic pressure, is
needed
(327)
The resulting gallls and timehistories for a variety of Mach numbers are described
further in Chapter 4
78

CHAPTER 4
RESULTS
It is the intent of the current effort to demonstrate the effectiveness of the
transpiration method in its application to steady and unsteady flow conditions. Based on
these results, the implementation of a discretetime control law within STARS is
discussed in regard to active fluttersuppression for the BACT wing. In a logical series of
steps, this chapter will present results starting with steadyflow simulations, which
include the effects of a deflected control surface, eventually leading up to both the open
and closedloop aeroservoelastic response. Where available, comparisons are made to
experimental data
4.] Steady Results Without Control Surface Deflections
The starting point of all unsteady cases in STARS, the steady state solution, must
be fully converged before starting an unsteady job For the final CFD mesh on the BACT
wing, the steady solution was mn for 3000 sieps to assure that the solution had, in fact,
converged to a steady state value. Convergence was assured using the maximum Mach
number criterion discussed in section 3.2.25. Experimental data are available for the
majority of test cases discussed in this section.
Remembering that the BACT CFD model is actually the CFD model for both the
NACA 0012 wing as well as the BACT wing, the differences between the two should be
79

noted here With an undeflected control surface, there should be no difference between
the two models since they are geometrically similar. Aside from slight manufacturing
differences, however, the two models were tested in a different fluid medium. The
NACA 0012 wing was tested in air (y=1.4) and the BACT wing was t.ested in R12
(y= I 148). As far as the calculation of the pressure coefficient is concerned, the value of
the ratio of specific heats, y, acts only as a scaling factor in steady simulations
Experimental steady data presented here comes from pressure transducers located at the
NACA 0012 wing's 60% span [Rivera, et ai, 1992]. More significant later, this 60%
span location corresponds to a distance of 19.2 inches from the wings root which, for the
BACT wing, corresponds to the midspan of the control surface.
The next six figures, Figure 41 to Figure 46, show steady pressure data obtained
at Mach numbers of 051, 0.67, 0.71, 077, 0.80, and 082, respectively. Each figure
shows pressure data at an angle of attack of 00
, with a control surface deflection of 0°
As each of the figures shows, agreement between predicted and experimental data is very
good, even at the higher transonic Mach numbers. Typically, as was briefly mentioned in
Chapter 2, Euler flow solvers overpredict both the location and strength of transonic
shocks. One common factor in each of the figures seems to be the fact that STARS tends
to predict a slightly higher suction peak, though still within the upper range of the
experimental scatter.
The BACT wing's critical Mach number appears to be 0.77 which coincides
with that of a NACA 00 12 airfoil At this point, flow accelerates from the freestream
Mach number of 0.77 to just sonic on the surface of the wing Cases run at Mach
numbers greater than 0 77 clearly show the existence of these transonic shocks
80

1.0
08
0.6
....... ~ ~.... .Ao ....
00.42 ..•.:.•. • • _~.., ......... . . ....... . .......
Goo • ~.
~&...
~2 • :
0.4 • Experiment
... STARS
0.6
0.8
0.5 0.6 0.7 08 0.9 1.0
xJc
0.1 0.2 0.3 04
1.0 .J.~~._r___~~I
00
Figure 41: Steady Chordwise Pressure at Mach 0.51, a.=0°, 8=0°, 60% Span
1.0
08 .
04
• Expenment
... STARS
0.6
0.8
09 10
1.0 +~~_r___~...._,r~t~~___+~~_'_t_~~_'__+~~~_+.~~.~r__~__+_~~_i
~o 01 02 OJ 04 05 ~6 07 08
xJc
Figure 42: Steady Chordwise Pressure at Mach 0.67, a=O°, 8=0°, 60% Span
81

1.0
·0.8
.(j 6
0.4
0.2
t •
8' 0.0 •
02
0.4 • Experiment
... STARS
0.6
0.8
1.0 +~~+'~~+~''>~I~~+~~'t~~~It""
~o ~I ~2 ~3 ~4 ~5 Q6 Q7 O~
x/e
0.9 1.0
Figure 43: Steady Chordwise Pressure at Mach 0.71, 0.=0°, 6=0°, 60% Span
1.0
0.8
0.4
06
• Experiment
... STARS
08
0.5 0.6 0.7 o.~ 09 J a
x/e
0.\ 02 0.3 04
I 0 l~~+'~~t'~''f'~I~~+~'t~~'+~~+~+~~1
00
Figure 44: Steady Chordwise Pressure at Mach 0.77, 0.=0°,6=0°,60% Span
82

1.0
04
•
• Experimenl
.A STARS
0.6
08
10
, I
08 09
'+'' 'f~" I
OA OS 06 07
xlc
0.1 0.2 0.3
1.0 +C~~+'~~t'~~t'~~t'~~
00
Figure 45: Steady Chordwise Pressure at Mach 0.80, 0.=0°, 6=0°, 60% Span
• Experiment
1.0
08
.06
04
.02
Co
U O.fl
.A STARS
OJ 04 05
x/e
06 07 08 O\! In
Figure 46. Steady Chordwise Pressure at Mach 0.82,0.=0°,6=0°,60% Span
...
83
b
Also accounting for the slight difference between predicted and experimental data
is a transition strip running approximately one inch from the leading edge of the wing.
There were no quantified uncertainties presented for these data, but judging from the
scatter in the pressure data, STARS predicts pressures that lie well within the
experimental scatter over the entire range of Mach numbers. Scatter is particularly
evident in Figure 45.
Solutions in the vicinity of Mach 0.77 took the most time to converge. At, and
slightly beyond, Mach 0.77 when shocks first begin to appear, solution convergence is
hampered as STARS resolves the location of the transonic shock. For lack of a better
term, the location of the shock seems to dance around a narrow portion of the wing's
surface. Though not a problem for the steady case, perse, a lack of resolution in the
shock locations could pose a problem with the unsteady flow solution. Addressed later,
the solution to this obstacle is to make sure that plenty of iterations are allowed for the
solution to completely converge at each solution timestep
4.2 Steady Results With Control Surface Deflections
The steady results presented above did not have to make use of the transpiration
boundary condition. For the case of a steady control surface deflection, there will be the
first actual application of the transpiration method thus far in this study. To show the
effectiveness of the transpiration method, a couple of different comparisons must be
made independent of one another. First, pressure distributions and contours are
compared for the case of a physical and transpired control surface deflection Second,
resulting pressure data for a simulated control surface deflection is compared to
experimental data horn the BACT wing
84
 
4.2.1 Steady Solutions for Transpired and Actual Control Surface Deflections
Recall that a CFD model for an actual control surface deflection was constructed
in addition to the standard CFD mesh for the wing. For purposes of comparison, a 100
control surface deflection is compared to that of a simulated 100 deflection angle. The
100 deflection angle was used to illustrate the effectiveness of tbe transpiration method
for relatively large surface deflections. Shown below in Figure 47 is a comparison
between the deflected and undeflected CFD grids.
'!I:1Aliilh'ImVsm:
l
p!l1 .m:eIlllICl!j!!l1 DDIiilU'ImWim·,!l1 .:i#!I§JOIIEISC1j
I
Figure 47: Comparison of Actual and Simulated 100 Control Surface Deflections
In the above figure, one can clearly see the extent of the flap deflection An Euler solver
would not be expected to detect or account for the likely separation and boundary layershock
interaction for such a large control surface deflection. The comparison with this
large control surface deflection is, therefore, used to demonstrate that the transpiration
method is as accurate as the limitations imposed by the inviscid flow assumption.
The first comparison of an actual and simulated control surface deflection is at
Mach 0.77, 00 angle of attack, and 100 (dowJJward) flap deflection A qualitative
comparison of the pictures in Figure 48 shows very good agreement between an actual
and simulated control surface deflection A more quantitative comparison can be made
85
with a comparison of the steady pressure distributions at the 60% span location, which
corresponds to the 'l2 span of the control surface.
10° ActliJial Deflection 100 Simulated Deflection
Figure 48: Surface Pressure Contours at Mach 0.77, 10° Control Surface Deflection
1 )
..... '..... I, .,~. . • Actual 100 Deflection • .._. r.....
0.4 • Simulaled 100 Deflection
0.8
0.6
,....... 
I ~ . I • •. •• ~ .4 I
"
I. .._ f.
0.4 . 4K I • '''.AI' I..... II 0.2 .. t .... ~
G °1
0.2 ..1
0.6 ..
05
xlc
0.8
I L._.~t'~~'~~~~"""'~ ,t'
o O. J 0.2 0.3 0.4 0.6 0.7 0.8 0.9
Figure 49: Comparison of Predicted Pressure Distributions for an Actual and Simulated
10° Control Surface Deflection at Mach 0.77, 0° (X
Figure 49 shows the excellent quantitative agreement between the predicted
pressure distribution for the actual and simulated control. surface deflection. With only
..
86
the slight discrepancy located at the x/c location which corresponds to the wing/control
surface interface. The rest of the data points essentially lie directly on top of oneanother.
The resulting differences in lift and moment predictions will also be small enough to be
considered insignificant.
At a slightly higher Mach number, Mach 0.82, similar results are presented. From
Figure 410 we again good qualitative agreement is seen between the pressure contours
not only on the wing, but out to the wall as well. Except for the fact that one can actually
see a physical deflection in the picture on the left, there is essentially no visual difference.
Quantitative agreement is again evaluated with the comparison of pressure distributions
at the 60% span location, see Figure 411. As was seen at Mach 0.77. the only noticeable
discrepancy between the pressure distributions is again at the same location, the
wing/control surface interface. One also notes the significant three dimensional effects
thal are captured as well.
: IiLplol VeI.ion 1 Db I!II!J E3
10° Actual Deflection
: IiLplol Vo..i"" 1.Db I!lIiIEJ
tJelp
10° Simulated Deflection

Figure 4] 0: Pressure Contours at Mach 0.82, 10° Control Surface Deflection
87
12
)
•...
0.8
.... _eJII ...... .. ~.. ~, ... £A.l 0.6 t ,. I 
0.4 f I \'.t :
021 • •...
G 0c,l '.~.' ~.~,
02 t .. Actual 10° Deflection I a. ........ a. .6t~
0.4 . • Simulated !00 Deflection • .: ..
0.6
0.8
0.1
t
I +1~~ro
.t02
0.3 0.4 os
x/c
06 07 0.8 09
Figure 411: Comparison of Predicted Pressure Distributions for an Actual and
Simulated 10° Control Surface Deflection at Mach 0.82, 0° a
The point as been made that the pressure distributions match well across the
chord, but what about the significant threedimensionality of the flow at the trailing edge.
With a control surface deflection, one would