February 2006 

NASA/TM-2006-214277 

Application of the ASP3D Computer Program 

to Unsteady Aerodynamic and Aeroelastic 

Analyses 

John T. Batina 

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

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role. 

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for NASA’s 

scientific and technical information. The NASA STI 

Program Office provides access to the NASA STI 

Database, the largest collection of aeronautical and 

space science STI in the world. The Program Office is 

also NASA’s institutional mechanism for 

disseminating the results of its research and 

development activities. These results are published by 

NASA in the NASA STI Report Series, which 

includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results ... even 

providing videos. 

For more information about the NASA STI Program 

Office, see the following: 

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 

at (301) 621-0134 

• Phone the NASA STI Help Desk at 

(301) 621-0390 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7121 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

February 2006 

NASA/TM-2006-214277 

Application of the ASP3D Computer Program 

to Unsteady Aerodynamic and Aeroelastic 

Analyses 

John T. Batina 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS) 

7121 Standard Drive 5285 Port Royal Road 

Hanover, MD 21076-1320 Springfield, VA 22161-2171 

(301) 621-0390 (703) 605-6000 
Application of the ASP3D Computer Program to 

Unsteady Aerodynamic and Aeroelastic Analyses 

John T. Batinaφ

NASA Langley Research Center 

Hampton, Virginia 23681 

Abstract 

A new computer program has been developed called ASP3D (Advanced Small 

Perturbation – 3D), which solves the small perturbation potential flow equation in an 

advanced form including mass-consistent surface and trailing wake boundary conditions, 

and entropy, vorticity, and viscous effects. The purpose of the program is for unsteady 

aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 

The program exploits the simplicity of stationary Cartesian meshes with the movement or 

deformation of the configuration under consideration incorporated into the solution 

algorithm through a planar surface boundary condition. The ASP3D code is the result of 

a decade of developmental work on improvements to the small perturbation formulation, 

performed while the author was employed as a Senior Research Scientist in the 

Configuration Aerodynamics Branch at the NASA Langley Research Center. The 

ASP3D code is a significant improvement to the state-of-the-art for transonic aeroelastic 

analyses over the CAP-TSD code (Computational Aeroelasticity Program – Transonic 

Small Disturbance), which was developed principally by the author in the mid-1980s. 

The paper presents unsteady aerodynamic and aeroelastic applications of ASP3D to 

assess the time dependent capability and demonstrate various features of the code. The 

cases considered include: (1) the NACA 0012 airfoil undergoing a forced harmonic 

pitching motion about the quarter chord at transonic conditions, (2) a thickening-thinning 

parabolic arc airfoil at a transonic Mach number, (3) an investigation of wave 

propagation characteristics with emphasis on high wave number applications, (4) 

aeroelastic transient calculations for the NACA 0012 airfoil at values of the dynamic 

pressure below, near, and above the flutter value, and (5) the F-5 fighter wing undergoing 

a rigid pitching motion about the wing root midchord axis. The results compare well 

with alternative methods and experimental data, and thus demonstrate the efficiency and 

utility of the ASP3D code for various unsteady aerodynamic and aeroelastic applications. 

Introduction

Over the years research has been conducted by several organizations including the 

NASA Langley Research Center to assess the general applicability and accuracy of 

transonic small perturbation (TSP) computer programs in numerous unsteady 

aerodynamic and aeroelastic applications.1-4 A notable example of such a computer 

 

φ

 Retired Research Consultant; Associate Fellow, American Institute of Aeronautics and 

Astronautics; formerly employed as a Senior Research Scientist, Configuration 

Aerodynamics Branch, NASA Langley Research Center, Hampton, Virginia. 
2

program is the CAP-TSD code5 developed by a team of researchers at NASA Langley 

Research Center in the mid-1980s. The CAP-TSD code solves the transonic small 

disturbance (TSD) equation with a time-accurate approximate factorization (AF1) 

algorithm6,7 and couples the flow equations with the structural equations of motion for 

simultaneous time integration for aeroelastic analyses. The code exploits the simplicity 

of stationary Cartesian meshes with the movement or deformation of the configuration 

under consideration incorporated into the solution algorithm through a planar surface 

boundary condition. However, applications with the code have revealed various 

inaccuracies inherent in the underlying small perturbation theory8

 including accuracy 

limitations near the leading edge and wing tip regions,1

 shock waves that are inaccurately 

captured in terms of strength or location,4

 and convergence difficulties attributable to 

certain computational and mathematical complications such as non-uniqueness of the 

potential flow equation.9,10 Consequently, the author conducted research over the last 

decade to examine the applicability of the small perturbation concept in general and the 

accuracy of the CAP-TSD code in specific. The author is the principal developer6,7,11,12 

of the CAP-TSD code and is therefore in a unique position to make a critical assessment 

of the methodology contained therein. The objective of the effort was to identify the 

sources of the inaccuracies and determine if improvements could be made to alleviate or 

eliminate them. 

Subsequently, a new advanced small perturbation (ASP) potential flow theory13 was 

developed by methodically determining the essential elements required to produce 

accurate solutions with the small perturbation approach on a Cartesian grid. The ASP 

theory involves a higher-order streamwise flux in the governing equation, massconsistent surface boundary conditions for flow tangency, mass-conserving entropy and 

vorticity effects including second-order terms in the trailing wake boundary condition, 

and viscous effects modeled using integral boundary layer equations14,15 solved 

simultaneously with the governing small perturbation potential flow equation with 

turbulent closure given by the dissipation integral relations of Drela.16 The ASP theory 

was shown to be mathematically more appropriate and computationally more accurate 

than the classical TSP theories.13 Therefore, the ASP theory was used as the basis for a 

new computer program called ASP3D17 (Advanced Small Perturbation – 3D), which 

involves either AF1- or AF2-type approximate factorization algorithms and a FAS (full 

approximation scheme) multigrid procedure for the solution of the ASP potential flow 

equation and associated boundary conditions. The purpose of the program is for unsteady 

aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 

The ASP3D code can treat aircraft configurations involving multiple lifting surfaces with 

leading and trailing edge control surfaces and a fuselage. Steady-state aerodynamic 

results obtained using the ASP3D code were presented by Batina,17 thus demonstrating 

various options and the general utility of the new program. Comparisons were made with 

alternative calculations and experimental data to demonstrate the accuracy of the ASP 

approach for steady-state applications. The purpose of the present paper, though, is to 

present unsteady aerodynamic and aeroelastic applications obtained using the ASP3D 

code. Several cases were considered to assess the time-dependent capability and 

demonstrate additional features of ASP3D. Comparisons are made with alternative 

calculations and experimental data to evaluate the accuracy and utility of the code. 
3

Overview of the ASP3D Program 

The ASP3D17 code was developed over the last decade by the author during his 

employment as a Senior Research Scientist within the Configuration Aerodynamics 

Branch at NASA Langley Research Center. The development of the code resulted from a 

research study to determine the specific causes of the inaccuracies and inefficiencies 

inherent to CAP-TSD. That work revealed that the underlying small perturbation 

potential flow theory in CAP-TSD was deficient in several respects. An advanced small 

perturbation theory13 was developed to alleviate or eliminate the inaccuracies identified 

in CAP-TSD. The ASP theory was shown subsequently to be mathematically more 

appropriate and computationally more accurate than the classical TSP theories such as 

that embodied within CAP-TSD. Consequently, the ASP theory was used as the basis for 

the new ASP3D program, which in contrast with CAP-TSD involves modern 

computational fluid dynamics concepts such as a finite volume spatial discretization, 

sharper shock capture, dissipation integral boundary layer modeling, alternative AF2-type 

approximate factorization algorithm, and FAS (full approximation scheme) multigrid 

solution convergence acceleration. 

The ASP3D Computer Program 

The ASP3D computer program involves either AF1- or AF2-type approximate 

factorization algorithms and a FAS multigrid procedure for the solution of the advance 

small perturbation potential flow equation and associated boundary conditions. The 

ASP3D code can treat aircraft configurations involving multiple lifting surfaces with 

leading and trailing edge control surfaces, and a fuselage. The new code is described 

briefly here along with relevant equations. The interested reader is referred to References 

13 and 17 for further details of the program. 

- Governing Equation 

The three-dimensional general small perturbation equation may be written in 

Cartesian coordinates as

∂f0

∂t

+

∂f1

∂x

+

∂f2

∂y

+

∂f3

∂z = 0

where the Cartesian fluxes are defined as 

f0 = −Aφ t − Bφ x

f1 = C + Dφx + Eφx

2 + Fφx3 + Gφy2 +

H

2 φxφ y

2

f2 = 1+ Hφx +

H

2 φx

2 ⎛

⎝

⎜ ⎞

⎠

⎟φ y + Fφy

3

f3 = φz 
4

with the constants defined as 

A = M∞

2 B = 2M∞2 C =1 D =1− M∞2

E = − 1

2

( ) γ +1 M∞

2 F = − 1

6

( ) γ +1 M∞

2 G = 1

2

( ) γ − 3 M∞

2 H = −( ) γ −1 M∞

2

The ASP3D code though solves the governing equation written in computational 

coordinates as 

∂ˆ f 0

∂t

+

∂ˆ f 1

∂ξ+

∂ˆ f 2

∂η

+

∂ˆ f 3

∂ζ = 0 

where the computational fluxes are defined as 

ˆ f 0 = − A

J

φt − B

J

ξ xφξ

ˆ f 1 = ∇ξ

J

ξ x

∇ξ f ( )1 + ξ y∇ξ f ( )2

⎡

⎣

⎢

⎤

⎦

⎥

ˆ f 2 = ∇η

J

ηy

∇η

f ( )2

⎡

⎣

⎢

⎤

⎦

⎥

ˆ f 3 = ∇ζ

J

ζ z

∇ζ f ( )3

⎡

⎣

⎢

⎤

⎦

⎥

The various geometric quantities used in these equations are computed using the exact 

formulas given by Batina.17

- Entropy and Vorticity Effects 

Shock-generated entropy effects are incorporated within small perturbation codes 

such as CAP-TSD by first using the Prandtl relation18 (shock jump condition) to 

determine the velocity downstream of the shock wave from the upstream and sonic 

velocities. The upstream and downstream velocities are then used in the RankineHugoniot shock relation18 to determine the change in entropy across the shock wave.19 

The resulting change in entropy is subsequently used in a Clebsch formulation20 to 

determine the vorticity downstream of the shock wave.21,22 The vorticity modifies the 

calculation of the velocity field in the downstream region. This is a common procedure 

that has been used in various computer codes especially at the full-potential equation 

level. However, the approach as applied to the small perturbation equation, with any of 

the available definitions for the streamwise mass flux f1 including that of Williams,23 does 

not conserve mass. 
5

In contrast, the approach developed for the ASP3D program conserves mass by 

using a ratio of the streamwise flux evaluated using the upstream and downstream 

velocities. Specifically, the downstream perturbation velocity is first computed using the 

Prandtl relation (shock jump relation) given by 

( ) φ x 2 = 1+ ( ) φx [ ] sonic

2

1+ ( ) φx 1

−1

where the subscripts 1 and 2 represent stations that are upstream and downstream of the 

shock wave, respectively. The change in entropy across the shock is then computed 

using 

Δs = (γ −1) 1− f1[ ] ( ) φ x 1 f { 1[ ] ( ) φx 2 }

For steady flows the entropy is held constant along gridlines downstream of shock waves. 

For unsteady flows the entropy is convected downstream using 

∂

∂t

Δs +

∂

∂x

Δs = 0 

The fluxes downstream of the shock are subsequently modified according to 

f ( )1 nonisentropic = 1− Δs

(γ −1)

⎡

⎣

⎢

⎤

⎦

⎥ f ( )1 isentropic

which conserves mass across the shock wave by design. 

The ASP3D code uses a Clebsch formulation20 to compute the shock-generated 

vorticity. In brief, the streamwise velocities downstream of shocks are computed using 

( ) φ x rotational = ( ) φx irrotational − Δs

γ(γ −1)M∞

2

Also, the second-order terms in the trailing wake boundary condition were found to be 

not insignificant, and thus, they were incorporated into the ASP3D solution procedure.17 

- Viscous Effects 

The ASP3D viscous method involves solving the unsteady boundary layer equations 

simultaneously with the outer potential flow solution so that no interaction law coupling 

the inner and outer solutions is required. The combined solution procedure involves an 

implicit block tridiagonal inversion for all of the cells along the surface and trailing wake. 

Unlike the CAP-TSD viscous capability,24-26 the ASP capability uses exact formulas for 

edge quantities and exact boundary conditions along surfaces and wakes. Smoothing of 
6

edge quantities and limiters are not required for stability, and no arbitrary or free 

parameters are necessary to tune the procedure. 

Specifically in the ASP3D viscous capability, the exact edge quantities are used for 

the calculation of velocity ue, density ρe, temperature Te, local Mach number Me, and the 

coefficient of viscosity μe defined by 

ue =1+ φx

ρe = 1− γ −1

2

M∞

2 2φx + φ x2 ( ) ⎡

⎣

⎢

⎤

⎦

⎥

1

γ−1

Te = 1− γ −1

2

M∞

2 2φx + φx2 ( ) ⎡

⎣

⎢

⎤

⎦

⎥

Me = M∞ ( ) 1+ φx

1− γ −1

2

M∞

2 2φx + φx

2 ( ) ⎡

⎣

⎢

⎤

⎦

⎥

1 2

μe = ( ) Te

3 2 1+ ( ) S T∞

( ) Te + ( ) S T∞

where S is Sutherland’s coefficient and T∞ is the freestream temperature. 

For the ASP3D capability, the time-dependent integral boundary layer (IBL) and lag 

equations may be written as a system of equations in the form 

[ ] A ∂

∂t

θ

Hk

cτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ [ ] B ∂

∂x

θ

Hk

cτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ [ ] C

φ xx

φxt

φtt

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ { }S = { }0

where the independent variables are the momentum thickness θ, the incompressible shape 

parameter Hk, and the square root of the shear stress coefficient cτ.

For solution, the integral boundary layer and lag equations are first premultiplied by 

[A]

-1 to diagonalize the time term as 

∂

∂t

θ

Hk

cτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ [ ] D ∂

∂x

θ

Hk

cτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ [ ] E

φxx

φxt

φtt

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ { }F = { }0
7

The resulting IBL system is then linearized in a simple way and cast into the so-called 

delta-form for implicit solution as 

Δx

Δθ

ΔHk

Δcτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ Δt G[ ]

Δθ

ΔHk

Δcτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

= −ΔxΔt R{ }

where the viscous residual is defined as 

{ }R = [ ] D ∂

∂x

θ

Hk

cτ

1 2

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ [ ] E

φxx

φxt

φ tt

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

+ { }F

The system of equations is directly coupled with the ASP equation for the outer 

potential flow for simultaneous implicit solution using a block tridiagonal matrix 

inversion procedure. This is done only for the cells adjacent to the upper and lower sides 

of the lifting surfaces and their wakes where the boundary layer is located. All other cells 

do not involve the boundary layer, and hence, only the ASP equation is solved in those 

cells. Dissipation integral relations are used for turbulent closure, similar to those 

reported by Drela.16

- Surface Boundary Condition 

In ASP3D the complete mass-consistent surface boundary condition including 

entropy, vorticity, and viscous effects is given by 

φz = f

g ( ) bx −α + φyby ± 1g

fδ

∗ ( )x + bt + δt*

where 

f = C + D( ) φx rotational + E φx

2 ( )rotational

+ F φx

3 ( )rotational

−Aφy

2 +

H

2 φy

2

( ) φx rotational

and 

g =1+ H( ) φx rotational+

H

2 φx

2 ( )rotational

+

H

2 φy

2 ( )

The surface is represented by 

b(x, y,t) ± δ

∗

(x,t) − xα = 0 

where the ± is for the upper and lower surfaces, respectively, since the boundary layer 

displacement thickness δ

*

 is defined as a positive quantity on both surfaces. 
8

Note that the spanwise slopes by are included in the complete surface boundary 

condition for three-dimensional applications, which may be important for swept wings 

near the leading edge and also in the wing tip region. Furthermore, the streamwise and 

spanwise slopes are computed internal to ASP3D from the ordinates of the wing 

geometry by simple second-order accurate central finite differences, effectively the same 

as is done in advanced codes with body fitted meshes. This eliminates the differentiation 

of the spline of the wing geometry normally performed as a preprocessing step to 

determine the airfoil slopes. It thus makes the input to the ASP3D code simpler than 

most other small perturbation codes such as CAP-TSD. 

- Wake Boundary Condition 

Viscous effects are included in the wake modeling in ASP3D through the boundary 

condition 

Δφz = ± 1

g

fδ

∗ ( )x

+ δt

∗

across the wake where f and g were defined previously. The boundary condition is 

incorporated numerically using 

( ) φ z upper = φk − ( ) φk−1 + Γ

zk − zk−1

+

1

2

1

g

fδ

∗ ( )x + δt∗ ⎡

⎣

⎢

⎤

⎦

⎥

and 

( ) φ z lower = ( ) φk − Γ − φk−1

zk − zk−1

− 1

2

1

g

fδ

∗ ( )x

+ δt

∗ ⎡

⎣

⎢

⎤

⎦

⎥

- Shock Capturing 

The AF1 and AF2 algorithms of the ASP3D code use one of three approaches to 

model supersonic regions of the flow including: (1) Godunov,27 (2) Engquist-Osher28 (EO), and (3) Murman-Cole29 (M-C) type-dependent switches. The first two switches 

satisfy an entropy inequality by design. This means that their use precludes the 

possibility of computing nonphysical expansion shocks that can occur with the MurmanCole switch, especially in two-dimensional applications. With Godunov27 switching, 

shocks are captured sharply, usually with only one interior state. Therefore, the Godunov 

switch is preferable over the E-O and M-C switches, and it is the default switch within 

the ASP3D code. 

- Solution Algorithm 

The ASP3D code involves a modern cell-centered finite-volume discretization, 

which allows exact planform treatment including leading edges, wing tips, and control 

surface edges. The advanced small perturbation equation is solved numerically using 

either an AF1 approximate factorization algorithm or an AF2-type30 algorithm. The AF1 
9

algorithm within the ASP3D code is the cell-centered finite-volume version of the AF1 

algorithm developed for CAP-TSD, which is described by Batina.7

 However, the AF1 

algorithm is susceptible to the so-called moving shock instability30 if shock waves move 

more than one cell per time step. 

The AF2 algorithm within the ASP3D code is described in general by 

LζLηLξΔφ = −R( ) φ

The general form of the AF2 algorithm appears to be similar to that of the AF1 algorithm 

but there are two fundamental differences. First, the operators are applied in the reverse 

order (here, vertical and spanwise sweeps are performed first, with the streamwise sweep 

performed last) and a streamwise implicit temporal damping term is included on the right 

hand side of the equation to eliminate the occurrence of the moving shock instability. 

The AF2 algorithm is defined specifically by 

α 

J

⎛

⎝

⎜ ⎞

⎠

⎟ −δζg ˆ 3

δζ

⎡

⎣

⎢

⎤

⎦

⎥Δφ i = −ωR +

1

J

ξ x

Δτ

+ψ

3ξ x

2Δt

B ⎛

⎝

⎜ ⎞

⎠

⎟Δφ i−1

α 

J

⎛

⎝

⎜ ⎞

⎠

⎟ −δηg ˆ 2δη

⎡

⎣

⎢

⎤

⎦

⎥Δφ i = α 

J

⎛

⎝

⎜ ⎞

⎠

⎟Δφ i

 

α 

J

⎛

⎝

⎜ ⎞

⎠

⎟ − g ˆ ( )1 i+1 2

r 

δ ξ

⎡

⎣

⎢

⎤

⎦

⎥Δφi = α 

J

⎛

⎝

⎜ ⎞

⎠

⎟Δφ i

where the flux Jacobians are defined by Batina.17 Although not obvious, the temporal 

damping term is treated implicitly since it involves information from the second 

intermediate sweep but at the previous grid plane i - 1, and 

α = ξ x

Δτ

+ψ

2

Δt

2 A +

3ξ x

2Δt

B ⎛

⎝

⎜ ⎞

⎠

⎟

The parameter ψ is set equal to zero for steady applications and set equal to unity for 

unsteady applications. The parameter ω is an over-relaxation factor. The AF2 algorithm 

involves both physical (Δt) and computational (Δτ) time steps. The computational step is 

selected for rapid convergence and it is usually assigned a large value. Hence the 

solution can be advanced using either a constant computational time step or a constant 

CFL number (variable time step). The physical time step is used only for unsteady 

applications, and it is selected based on the physics of the problem being considered. 

The AF2 approach is the preferred algorithm. The AF2 scheme is more robust than 

the AF1 scheme since it cannot fail due to the moving shock instability by design, 

because the temporal damping term eliminates the instability. Hence convergence to 

steady state is achieved with large time steps or large CFL numbers, independent of the 
10

moving shock instability that plagues the AF1 algorithm. The AF2 scheme also may be 

used for unsteady calculations when local iteration is applied. Therefore, the physical 

time step may be selected based on the physics of the problem being solved, rather than 

on numerical stability considerations. In either steady or unsteady applications, the 

temporal damping term vanishes in the steady-state limit or within the convergence 

process of performing subiterations (for unsteady calculations). 

- Multigrid Implementation 

The AF1 algorithm is not appropriate for multigrid implementation because the rapid 

rate of convergence of the multigrid procedure would result in the moving shock 

instability. The AF2 algorithm of the ASP3D code is amenable to multigrid 

implementation because the temporal damping term prevents the moving shock 

instability from occurring. Therefore the multigrid procedure31 may be used for steady or 

iterative unsteady applications with ASP3D. The procedure is implemented in FAS (full 

approximation scheme) form,32 which is applicable to nonlinear governing equations. As 

with any multigrid procedure, high-frequency errors are damped on the fine mesh and 

low-frequency errors are damped on the coarser meshes. Either V- or W-cycles may be 

selected and full multigrid is available. W-cycles have been found to be slightly more 

effective for solution convergence, although W-cycles are slightly more expensive than 

V-cycles.33

Specifically, the calculation on the finest mesh (N) smoothes high frequency errors 

using the AF2 algorithm written in general form 

LN ( ) φ N = RN

Subsequent calculations on the next coarser mesh (N-1) involve fine mesh residual 

injection 

LN −1( ) φ N−1 = RN−1 − ˆ I N

N−1

( ) RN + LN−1IN

N−1

( ) φ N

where the I functions are restriction operators for transferring the velocity potential and 

residual from the fine mesh to the coarse mesh. The procedure is generally repeated for 

additional coarse meshes (N-2, N-3, ...). 

The restriction operator33 for the velocity potential is a volume-weighted operator 

defined by 

φ N−1 = IN

N−1

φ N = ∑Vφ N

∑V

and the restriction operator33 for residuals is defined by 

ˆ I N

N−1

RN = ∑ RN

where the summations above are taken over all of the fine mesh cells that make up the 

coarse mesh cell. The interpolation operator for the velocity potential that is required to 
11

transfer information back from coarse to fine mesh is a trilinear interpolation operator 

defined symbolically by 

φ N ← φ N + IN−1

N φ N−1

- Numerical Accuracy 

In the ASP3D code geometric quantities such as areas, volumes, and direction 

cosines are computed using exact formulas.17 The velocity components are computed 

using similar formulas as used in CAP-TSD, but in ASP3D, they are exact for a linear 

test function on a general grid. Generally speaking, the ASP3D algorithm is secondorder-accurate within subsonic cells and first-order accurate within supersonic cells, 

although a second-order-accurate supersonic capability is available as a user-selected 

option. 

Also in CAP-TSD the φxt term in the governing equation is computed using firstorder-accurate backward differencing for stability reasons. However, for unsteady 

applications at higher reduced frequencies encountered with control surfaces or higher 

wave number phenomena such as Kutta waves, this differencing may be too dispersive. 

In contrast, first-, second-, and third-order-accurate backward differencing operators are 

available for the φxt term in ASP3D as user-selected options to improve the accuracy for 

such applications. 

- Structural Equations of Motion 

Similar to the CAP-TSD code, the flow equations of ASP3D are coupled with the 

structural equations of motion for simultaneous time-integration for aeroelastic analysis. 

Specifically, such analysis involves the coupling of the aerodynamics with the structural 

characteristics of the configuration under consideration.2

 The resulting equations of 

motion for a time domain or time-marching aeroelastic solution are based on the aircraft 

natural vibration modes. These equations are integrated in time along with the finite 

volume solution of the flow field in ASP3D. Initial conditions for each structural mode 

are specified and free decay transients are computed. Initial velocities in one or more 

modes, rather than displacements, have been found to be superior in avoiding 

nonphysical, strictly numerical transients and their possibly associated instabilities. 

Aeroelastic stability is then deduced from the free decay records or time histories. This is 

a fairly standard procedure for aeroelastic analyses with the small perturbation computer 

codes, and thus, the interested reader is referred to Cunningham, et al.2

 or Bennett, et al.3 

for further details including equations and representative calculations. 

Results and Discussion 

Results computed using the ASP3D computer program are presented for several 

unsteady aerodynamic and aeroelastic cases to assess the time-dependent capability and 

demonstrate various features of the code. The cases considered include: (1) the NACA 

0012 airfoil undergoing a forced harmonic pitching motion about the quarter chord at 
12

transonic conditions,34 (2) a thickening-thinning parabolic arc airfoil at a transonic Mach 

number,35 (3) an investigation of wave propagation characteristics with emphasis on high 

wave number applications,36,37 (4) aeroelastic transient calculations for the NACA 0012 

airfoil at values of the dynamic pressure below, near, and above the flutter value,38,39 and 

(5) the F-5 fighter wing undergoing a rigid pitching motion about the wing root midchord axis.40 All of the calculations are inviscid applications of the code, in part, because 

the viscous capability was not implemented at the time. 

Flow calculations for the two-dimensional cases were performed using a 193 x 129 

Cartesian finite volume mesh as shown in Figure 1. The total mesh is shown in Figure 

1(a) and a near field view is shown in Figure 1(b). The mesh had a twenty chord extent 

from the airfoil in both the x- and z-directions with a spacing at the leading and trailing 

edges of Δx = 0.015c. There is only one cell in the spanwise direction for airfoil 

applications, which effectively produces a two-dimensional solution. Multigrid 

calculations were performed using this mesh (termed Mesh 5) and four subsets of the 

mesh (Meshes 1 through 4) determined by simply deleting every other grid line. The 

coarser meshes contain 97 x 65 (Mesh 4), 49 x 33 (Mesh 3), 25 x 17 (Mesh 2), and 13 x 9 

(Mesh 1) points in the x- and z-directions, respectively. This is done automatically 

within the ASP3D code and thus is not something that the user needs to do during mesh 

generation. Flow calculations for the F-5 fighter wing were performed using a 97 x 25 x 

65 finite volume mesh as discussed below. 

- NACA 0012 Airfoil Harmonic Pitching Results

The first case considered is that of the NACA 0012 airfoil at a mean angle of attack 

of α0 = 0.016° and a freestream Mach number of M∞ = 0.755.34 The airfoil was forced to 

pitch harmonically about the quarter chord with an amplitude of α1 = 2.51° at a reduced 

frequency based on semichord of k = 0.0814. This is a well-studied nonlinear unsteady 

transonic case34 that involves moving shock waves on both the upper and lower surfaces 

of the airfoil produced by the relatively large amplitude of motion. The shock waves 

move approximately twenty-five percent chord and periodically appear and disappear 

during a cycle of motion. The case was chosen to demonstrate the superior numerical 

stability of the iterative multigrid capability of the ASP3D code for unsteady 

applications, whereby the computational time step is selected to be very large for rapid 

convergence and the physical time step is selected based on the physics of the problem 

(temporal accuracy). In fact the physical step size can be selected to be very large as 

well, without creating numerical difficulties. 

To demonstrate the algorithm robustness, iterative multigrid calculations were 

performed for the NACA 0012 airfoil oscillating harmonically for four cycles of motion 

to ensure periodicity of the last cycle. The CFL number was selected to be 10 (which 

corresponds to a spatially varying computational step size) and the over-relaxation factor 

(ω) was set equal to 1.6. The convergence criteria on a per time step basis was selected 

as two orders of magnitude, which produced an L2-norm of the residual of approximately 

10-3 at each time step. Calculations were performed first using only 32 steps per cycle of 

motion corresponding to a relatively large physical step size of Δt = 1.20608. This step 
13

size corresponds to more than one chord of airfoil travel per time step. The forced 

pitching motion α(t) is shown in Figure 2(a) and the lift and moment coefficient 

responses, cl(t) and cm(t), respectively, are shown together in Figure 2(b). The 

coefficients indicate that the lift is a relatively linear response that slightly lags the 

pitching motion and the moment is a highly nonlinear response that lags the motion much 

more than the lift. Figures 2(c) and 2(d) show the iterative multigrid residual history and 

the number of supersonic points, respectively, during the last cycle of motion. The 

residual history essentially represents thirty-two successive multigrid calculations, one 

per physical time step, wherein the multigrid procedure is used to rapidly reduce the 

solution residual to the desired level, thus preserving temporal accuracy. By iterating on 

the solution at each time step, the linearization and factorization errors inherent in the 

AF2 algorithm are reduced to an acceptable level (~10-3) to produce a time accurate 

solution. The technique is analogous to performing subiterations per time step with CAPTSD to ensure time accuracy, but the iterative multigrid procedure in ASP3D is cheaper 

computationally because the iterations are performed on coarser grids within the overall 

multigrid capability. The ASP3D calculation required about 180 iterations of the 

multigrid procedure, but since most of the work is done on the coarser Meshes 1 through 

4 using W-cycles, rather than on the finest Mesh 5, the total effort is less than iterating 

using only the finest mesh. The number of supersonic points during the last cycle of 

pitching motion (Figure 2(d)) is cyclic, as expected, with a larger value in the first half of 

the cycle than the value attained in the second half because of the small positive mean 

angle of attack (α0 = 0.016°). 

To further test and demonstrate the algorithm robustness of the ASP3D code, 

calculations for the pitching NACA 0012 airfoil were successively repeated using 

multiples of the original physical step size including Δt = 2.41216, 4.82432, and 9.64864, 

which correspond to 16, 8, and 4 steps per cycle (s/c) of motion, respectively. These 

values represent huge step sizes that are not usually numerically stable with other TSP 

computer codes. Admittedly, 16 and 8 steps per cycle barely describe the pitching 

motion in discrete form, and 4 steps per cycle of motion is the smallest value possible, 

which actually corresponds to a saw tooth pitching motion rather than a sinusoidal 

motion. Be that as it may, all of the calculations were numerically stable, and the results 

for 16, 8, and 4 s/c are presented in Figures 3, 4, and 5. The results are presented in the 

same format as the results of Figure 2 computed using 32 steps per cycle. Specifically, 

the pitching motion, the lift and moment coefficient time responses, the last cycle 

multigrid residual history, and the last cycle number of supersonic points are given in 

sub-figures (a), (b), (c), and (d), respectively. Comparisons of the results from the 16 s/c 

calculation (Figure 3) with the results obtained using 32 s/c (Figure 2) indicate that the 

lift response is very similar and the moment response is not as smooth, possibly because 

it is difficult to plot the nonlinear moment time history with only 16 straight line 

segments per cycle. The last cycle residual history shows very similar convergence to the 

32 s/c case but only 96 iterations of the multigrid procedure were required to obtain 

results for 16 s/c, which is approximately half the work required for the 32 s/c results. 

Comparisons of the results obtained using 8 s/c (Figures 4(a) – 4(d)) and 4 s/c 

(Figures 5(a) – 5(d)) with those using 32 s/c (Figures 2(a) – 2(d)) indicate a degradation 
14

in the accuracy of the responses, especially the moment coefficient, but the calculations 

are numerically stable even though the corresponding step sizes are huge. The number of 

multigrid iterations required to obtain results with 8 s/c and 4 s/c are 48 and 24, 

respectively, as indicated in Figures 4(c) and 5(c). These values are exactly half and onequarter the work required to obtain the results with 16 s/c (96 iterations), which indicates 

that there is a linear relationship between the number of steps per cycle and the number of 

multigrid iterations required to obtain results. For example, halving the number of steps 

per cycle in the calculation cuts the total work in half. 

Composite hysteresis plots of the responses in the last cycle of motion are plotted in 

Figure 6. Figure 6(a) shows the lift coefficient cl(t) versus the pitching motion α(t) and 

Figure 6(b) shows the moment coefficient cm(t) versus the pitching motion α(t). The 

results from all four sets of calculations are presented corresponding to 32, 16, 8, and 4 

s/c of motion. The responses obtained using 32 and 16 s/c compare closely (the discrete 

points are connected with straight line segments) while the responses obtained using 8 

and 4 s/c are less accurate because the corresponding step sizes are too large. However, 

for easier cases involving smaller amplitudes of motion and less shock motion including 

flutter onset calculations (due to a small perturbation of the flow) the larger time steps 

may produce solutions of acceptable accuracy at a fraction of the cost. In any case, use 

of larger time steps does not create numerical stability problems, and therefore the analyst 

must make a prudent selection of step size balancing accuracy and cost of calculation. 

- Thickening-Thinning Parabolic Arc Airfoil Results

The second case considered is that of a parabolic arc airfoil with a time-varying 

midchord thickness at M∞ = 0.85 and α0 = 0°.

35 The temporal variation in thickness is 

defined mathematically and graphically in Figure 7. For example, beginning as a flat 

plate at t = 0, the thickness of the parabolic arc airfoil increases smoothly until it reaches 

maximum thickness of 10% (maximum thickness to chord ratio) at t = 15. After t = 15 

the airfoil thins until it reaches zero thickness at t = 30. This case is a simple twodimensional helicopter rotor simulation that has been used by other researchers for code 

development purposes.35 Similarly in the present study, the problem was solved using 

the ASP3D iterative multigrid procedure involving the five meshes discussed previously. 

The calculations were performed using a CFL number of 10 and a relatively large 

physical step size of Δt = 1.0 for a total of 32 steps. The step size of Δt = 1.0 corresponds 

to one chord of airfoil travel per time step. The convergence criteria on a per time step 

basis was selected as four orders of magnitude, which produced an L2-norm of the 

residual of approximately 10-5 at each time step. The resulting residual history and 

number of supersonic points are shown in Figures 8(a) and 8(b), respectively. The 

residual history (Figure 8(a)) indicates that 184 iterations were required to advance the 

solution 32 time steps and that each time step required about the same amount of work (5 

or 6 iterations) independent of whether the flow was subsonic or transonic. The number 

of supersonic points (Figure 8(b)) indicates that the solution goes from subsonic to 

transonic flow at time step 7, with a maximum of about 1250 supersonic points near time 

step 18 or 19, and back to subsonic flow at time step 27. 
15

For comparison purposes the calculations were repeated with everything the same 

except for using only the finest Mesh 5. This calculation is referred to as an iterative 

unigrid calculation (no multigrid) to contrast it with the iterative multigrid calculation of 

Figure 8. The procedure then only iterates with the solution on the finest mesh, similar to 

that which is done in the CAP-TSD code when subiterations are applied on a per time 

step basis. The resulting residual history and number of supersonic points are presented 

in Figures 9(a) and 9(b), respectively. The residual history (Figure 9(a)) indicates that 

over 2400 unigrid iterations were required to achieve a similar convergence level as the 

iterative multigrid calculation (L2-norm of the residual of approximately 10-5 at each time 

step) and clearly more work (number of iterations) is required when the flow is transonic 

than when the flow is subsonic. Thus the iterative multigrid calculation (184 total 

iterations) is far more efficient than the iterative unigrid calculation (over 2400 total 

iterations), because most of the iterative multigrid work is done on coarser meshes 

(Meshes 1 through 4). The number of supersonic points (Figure 9(b)) is similar to that 

shown from the multigrid calculation (Figure 8(b)), such as attaining a maximum of 

about 1250 supersonic points, since the two flow field solutions should be almost 

identical. 

Instantaneous solutions from the iterative multigrid calculations at t = 8.5, 11.5, 

18.25, 26.875, 29.0, and 32.0 are shown in Figures 10, 11, 12, 13, 14, and 15, 

respectively. These six points in time were selected for later comparison with the 

published results of Goorjian.35 The top part (a) of each figure shows the surface 

pressure coefficient distribution and the bottom part (b) of each figure shows the near 

field pressure coefficient contour lines. For t = 29.0 and 32.0, the results were taken 

directly from the multigrid calculations discussed above since the step size was Δt = 1.0 

(exactly 29 and 32 total time steps, respectively). To avoid solution interpolation to 

obtain the results for t = 8.5, 11.5, 18.25, and 26.875, the calculations were repeated four 

times with step sizes equal to Δt = 8.5/8, 11.5/11, 18.25/18, and 26.875/26, respectively, 

using 8, 11, 18, and 26 total time steps. 

At t = 0 the airfoil has zero thickness (flat plate) so the undisturbed flow is the 

correct starting solution for the thickening-thinning airfoil problem. As the thickness of 

the parabolic arc airfoil increases in time the solution develops such that, for example at t

= 8.5 (Figures 10(a) and 10(b)), the flow is relatively smooth with very small embedded 

supersonic regions forming on each side of the airfoil with weak shock waves beginning 

to develop near 76% chord. As the airfoil further thickens to t = 11.5 (Figures 11(a) and 

11(b)) the shock waves grow to moderate strength near 83% chord. The airfoil reaches 

its maximum thickness of 10% at t = 15, so by t = 18.25 (Figures 12(a) and 12(b)), the 

embedded supersonic regions are approximately at their largest (1250 supersonic points) 

and the relatively strong shock waves are located near 91% chord. As the airfoil becomes 

thinner in time the flow becomes subsonic again and the upper and lower surface shock 

waves become upstream propagating pressure waves near midchord as indicated at t = 

26.875 (Figures 13(a) and 13(b)). At t = 29.0 (Figures 14(a) and 14(b)) the pressure 

waves have traveled upstream to roughly the airfoil quarter-chord location and by t = 

32.0 (Figures 15(a) and 15(b)), when the airfoil has zero thickness again, the pressure 
16

waves on the upper and lower surfaces have moved forward of the airfoil leading edge 

thus forming a single wave. 

Comparisons of surface pressure coefficients for the thickening-thinning parabolic 

arc airfoil case are presented in Figure 16. The ASP3D pressure distributions at t = 8.5, 

11.5, and 18.25 are shown in Figure 16(a) with corresponding full-potential (FP) and 

TSD pressure distributions from Goorjian35 shown in Figure 16(b). The FP and TSD 

results of Reference 35 were obtained using step sizes of 0.0625 and 0.125, respectively, 

whereas the ASP3D results were obtained using step sizes equal to or greater than 1.0. 

The ASP3D pressures (Figure 16(a)) are in good agreement with the pressures reported 

by Goorjjian,35 (Figure 16(b)) especially in comparison with the full-potential pressure 

results. Specifically at t = 18.25 where the shock waves are strongest, the ASP3D shock 

location compares better with the FP result than the TSD result. Further pressure 

comparisons for t = 26.875, 29.0, and 32.0 are presented in Figure 16(c) from ASP3D 

and in Figure 16(d) from Reference 35. The ASP3D pressures again agree best with the 

full-potential pressure distributions, although all of the results indicate that the pressure 

waves are smeared as they propagate upstream off the airfoil leading edge. 

- Linear Wave Propagation Results 

To understand the importance of accurate wave propagation for unsteady 

aerodynamic flows, it is instructive to relate some of the work reported over thirty years 

ago by Tijdeman and Zwaan36 at the NLR in Amsterdam. Reference 36 contains a 

discussion of the propagation of pressure waves in high subsonic and transonic flows. 

The authors looked at disturbance wave fronts generated by a harmonically pulsating 

source at the three-quarter-chord location (control surface hinge line) along an airfoil, 

presented here in Figure 17 taken from their report. Figure 17(a) shows upstream 

propagating wave fronts for a subsonic freestream Mach number of M∞ = 0.8 and Figure 

17(b) similarly shows wave fronts for a transonic Mach number of M∞ = 0.875. The 

latter flow involves an embedded supersonic region on the airfoil terminated by a shock 

wave near midchord. In the upper parts of the two figures the positions of the wave 

fronts are drawn at equal time intervals for the slowly moving receding waves. The 

lower parts of the figures show the same wave fronts, but for a uniform flow (i.e., without 

the airfoil) as a reference. At the subsonic freestream Mach number of M∞ = 0.8 (Figure 

17(a)), the actual (upper) wave fronts move a little slower than the uniform flow (lower) 

as suggested by the closer spacing of the curves representing the fronts. The spacing of 

the wave fronts is a measure of the intensity of the local pressure gradients. Also the 

wave fronts incline forward slightly due to the vertical velocity gradient. At the transonic 

Mach number of M∞ = 0.875 (Figure 17(b)), the lower part of the actual wave fronts 

coalesce into the shock wave. The upper parts of the wave fronts rotate about the top of 

the shock and enter the supersonic region at an oblique angle, inclined in comparison to 

the uniform flow. The waves are slowed significantly near the shock wave thus 

producing a strong pressure gradient there, as expected. Also, the longer time required 

for the wave fronts to travel around the shock wave leads to a significant time lag in the 

airfoil aerodynamic response. Although not demonstrated here, this time lag is generally 

regarded as the physical mechanism responsible for the so-called transonic flutter dip. 
17

Tijdeman and Zwaan36 also presented wind tunnel results in the form of 

shadowgraph photographs of high-speed flows about an airfoil reproduced here in Figure 

18. The upstream propagating wave fronts are generated by vortex streets emanating 

from the aft end of the airfoil with its trailing edge control surface removed. Figure 18(a) 

shows the shadowgraph for the subsonic Mach number of M∞ = 0.8 wherein the receding 

wave fronts have relatively uniform spacing and they are slightly inclined as in the 

diagram of Figure 17(a). Figure 18(b) shows the shadowgraph for the transonic Mach 

number of M∞ = 0.875 wherein the lower parts of the receding wave fronts coalesce near 

the foot of the shock wave and the upper parts of the fronts rotate around the top of the 

shock to enter the supersonic region, similar to that diagrammed in Figure 17(b). These 

figures illustrate some of the flow field complexities for unsteady transonic aerodynamics 

that computational methods must be able to model accurately. 

It also is instructive to examine the simpler linear theory wave propagation 

characteristics of the upstream and downstream propagating wave fronts in a uniform 

flow and the implications for resolving such waves using discrete computational 

methods. For example, Figures 19(a) and 19(b) show wave numbers for upstream and 

downstream propagating wave fronts, respectively, as a function of Mach number for 

various values of reduced frequency k. The wave number associated with upstream 

propagating wave fronts is defined by 

 wave number = M∞

1− M∞

k

and the wave number associated with downstream propagating wave fronts is defined by 

 wave number = M∞

1+ M∞

k

According to Rowe and Cunningham,37 “The wave number is a useful parameter for 

assessing the waviness of pressure distributions and estimating relative difficulties in 

predicting converged values of generalized forces.” Computational fluid dynamics 

(CFD) methods such as CAP-TSD or ASP3D are typically capable of computing 

propagating wave fronts for wave numbers less than about 5.0 depending upon the 

density of the mesh and the size of the physical time step. Therefore the computation of 

most downstream propagating waves (shown in Figure 19(b)) with associated wave 

numbers less than about 2.0 is not difficult because the waves are moving very fast (with 

the freestream), although relatively small time steps may be required to accurately 

describe downstream running waves at higher reduced frequency. In contrast, the 

computation of upstream propagating wave fronts is difficult at higher reduced 

frequencies because the associated wave number is much higher, as shown in Figure 

19(a). Pressure distributions with wave numbers higher than about 5.0 have significant 

waviness that can be difficult to capture with CFD methods as discussed below. As the 

Mach number increases toward unity from the subsonic side the problem is exacerbated, 

because the waves are moving increasingly slowly (against the freestream). The wave 

number curves in Figure 19(a) asymptote to infinity at Mach one because the wave fronts 
18

there are theoretically standing waves. And for supersonic freestream there can be no 

upstream propagation. 

Another way of looking at the linear theory wave propagation characteristics is by 

plotting the reduced frequency of the upstream and downstream propagating wave fronts 

as functions of Mach number for various values of the wavelength λ as presented in 

Figures 20(a) and 20(b), respectively. The upstream values are computed by 

k = 1− M∞

M∞

π

λ

 

and the downstream values are computed by 

 k = 1+ M∞

M∞

π

λ

 

For the upstream propagating wave fronts (Figure 20(a)) the reduced frequency and 

wavelength are inversely proportional. As the Mach number approaches unity from the 

subsonic side the reduced frequency is reduced for a given wavelength until the singular 

limit at Mach one is reached where a standing wave can occur. The small wavelength 

fronts that may be regarded as being outside the range of practical reduced frequencies 

for Mach numbers of, say, M∞ = 0.8, are in fact within the range encountered in 

aeroelasticity (approximately 0 < k < 2) at higher Mach numbers of 0.9 or 0.95. 

Therefore the accurate resolution of these small wavelength wave fronts should be of 

concern for CFD development. Of course for supersonic freestream there is no upstream 

wave propagation. For the downstream propagating wave fronts (Figure 20(b)) the 

reduced frequency is again inversely proportional to the wavelength. But here, 

wavelengths less than about 4.0 result in reduced frequencies that are outside the typical 

range encountered in aeroelasticity. Conversely, wavelengths greater than 4.0 correspond 

to reduced frequencies less than 2.0 for any Mach number, but since these wavelengths 

are so large, they can be represented easily with most CFD codes given a sufficient grid 

and appropriate step size. Therefore it is the upstream propagating waves where the 

performance of these codes needs to be tested.

Rowe and Cunningham37 presented lifting pressure coefficient distributions for 

unsteady cases involving high wave numbers. Similarly, Figure 21 shows the effects of 

Mach number on the lifting pressure coefficient distributions for a flat plate airfoil 

oscillating in pitch about the leading edge at k = 0.9. Figure 21(a) shows the real or inphase component of the pressure and Figure 21(b) shows the imaginary or out-of-phase 

component of the pressure. The results were computed using the two-dimensional 

assumed-pressure-mode kernel-function linear-aerodynamics program of Bland.41 The 

figures clearly show that as the freestream Mach number increases from 0.75 to 0.95 

more waves occur in the pressure distribution. For the Mach numbers of 0.75, 0.80, 0.85, 

0.90, and 0.95, the corresponding upstream propagating wave numbers are 2.7, 3.6, 5.1, 

8.1, and 19.0, respectively. 
19

Figure 22(a) shows the lifting pressure coefficient distribution computed using the 

ASP3D code for the flat plate airfoil oscillating in pitch about the leading edge at k = 0.9. 

The freestream Mach number was selected as M∞ = 0.9 and the linear theory results from 

the previous Figure 21 that had a relatively high wave number of 8.1 are repeated here for 

comparison purposes. It also is a case studied in Reference 37. Figure 22(a) shows the 

real and imaginary ASP3D pressures, which agree with the linear theory pressures near 

the trailing edge, but do not agree with linear theory as the waves propagate upstream 

toward the leading edge. The ASP3D results were obtained using a first-order-accurate 

backward difference formula to compute the φxt term. A backward difference formula is 

typically used to compute the φxt term in TSP codes for stability reasons, and similar 

erroneous results were obtained using CAP-TSD. However in ASP3D, second- and 

third-order-accurate formulas have been successfully implemented to improve the 

temporal accuracy of the results for higher wave number cases. The pressure distribution 

for this case computed using the second-order-accurate formula is presented in Figure 

22(b). These results show a marked improvement in accuracy in comparison with the 

linear theory result. Results using the third-order-accurate formula are identical to 

plotting accuracy with the second-order-accurate results, and thus, they are not presented 

here. 

- NACA 0012 Airfoil Aeroelastic Transient Results 

To assess and demonstrate the aeroelasticity capability in ASP3D, calculations were 

performed for the NACA 0012 airfoil at M∞ = 0.8 and α0 = 0°, a case considered by 

Robinson, et al.38 using the CFL3D42 code run in Euler mode. Steady-state calculations 

were performed first, using ASP3D including entropy and vorticity effects, to provide the 

initial conditions for aeroelastic computation upon restart. The resulting steady pressure 

coefficient distribution is presented in Figure 23(a) with the corresponding CFL3D 

(Euler) pressure distribution from Ref. 38 presented in Figure 23(b). The ASP3D 

pressure distribution is in very good agreement with the CFL3D pressure distribution 

whereby the shock waves on the upper and lower surfaces of the airfoil are predicted to 

be near midchord. In the ASP3D calculation the shock waves are captured sharply 

because of the use of the Godunov type-dependent switch and the pressure levels are well 

predicted including those near the airfoil leading edge (stagnation). 

For aeroelastic analysis, the NACA 0012 airfoil was modeled structurally with two 

degrees of freedom, bending and torsion. The case considered is the much-studied Case 

“A” of Isogai39 which has normal modes similar to those of a streamwise section of a 

sweptback wing. The wind-off bending and torsion natural frequencies are 71.33 and 

535.65 rad/s, respectively. The pivot point for the bending mode is located 1.44 

chordlengths upstream of the leading edge of the airfoil. The pivot point for the torsion 

mode is 0.068 chordlengths forward of midchord. These mode shapes and natural 

frequencies were determined by performing a free vibration analysis with the aeroelastic 

equations written in the traditional form of plunge and pitch degrees of freedom. In this 

analysis, the following structural parameter values were used: a = -2.0, xα = 1.8, rα = 

1.865, ωh = 100 rad/s, and ωα = 100 rad/s. Also, the airfoil mass ratio was μ = 60. 

Generalized displacements corresponding to the bending and torsion modes are defined 
20

as q1(t) and q2(t), respectively. Initial conditions for the time-marching aeroelastic 

calculations were dq1(0)/dt = 2.0 and dq2(0)/dt = 1.0. 

Aeroelastic calculations were performed for the NACA 0012 airfoil at M∞ = 0.8 and 

α0 = 0° for several values of the nondimensional dynamic pressure including Q = 0.2, 

0.5, and 0.8, where Q = (U∞/bωα√μ)

2

. In these calculations a smaller value of the step 

size was selected than in previous unsteady calculations, Δt = 0.1, to ensure that the 

higher frequency torsion mode would be resolved in the aeroelastic integration. Also the 

multigrid procedure was not exercised because the computational benefit would not be as 

large for such a step size. The generalized displacements for Q = 0.2, 0.5, and 0.8 are 

presented in Figures 24, 25, and 26, respectively. For Q = 0.2 (Figure 24), the 

generalized displacements indicate that the aeroelastic system is stable because the 

transients are converging in time. The q1(t) generalized displacement is dominated by the 

lower frequency bending mode and the q2(t) generalized displacement contains both 

modes including the higher frequency torsion mode. Use of the step size of Δt = 0.1 

resulted in approximately 133 s/c for the bending mode and only 20 s/c for the higher 

frequency torsion mode. Use of much larger time steps, although numerically stable, 

produces results that inaccurately describe the torsion mode. 

For Q = 0.5 (Figure 25), the generalized displacements indicate that the aeroelastic 

system is near-neutrally stable (flutter) since the dominant bending mode in both q1(t)

and q2(t) is not damped. This correlates well with the reported results from Ref. 38 

wherein the transients computed using CFL3D (Euler) at Q = 0.5 were described as nearneutrally stable corresponding to flutter. Also the q2(t) generalized displacement 

indicates that the higher frequency torsion mode is more damped than at Q = 0.2. Finally 

for Q = 0.8 (Figure 26), the generalized displacements q1(t) and q2(t) are both diverging 

corresponding to an unstable aeroelastic system and the higher frequency torsion mode 

contained within q2(t) has even more damping than the Q = 0.2 and 0.5 transients. 

- F-5 Fighter Wing Harmonic Pitching Results

The F-5 fighter wing40 has a leading edge sweep of 31.92 degrees, an aspect ratio of 

1.578, and a taper ratio of 0.283. The wing has a NACA 65A004.8 airfoil section that 

has a drooped leading edge. The wing was modeled for ASP3D applications using a 97 x 

25 x 65 finite volume mesh as shown in Figure 27. A near field view of the 97 x 25 

planform mesh is shown in Figure 27(a) and a near field view of the 97 x 65 root 

sectional mesh is shown in Figure 27(b). The meshes were generated using a program43

that was created to construct meshes specifically for use by the ASP3D code, and as such, 

it can generate meshes for aircraft configurations with multiple lifting surfaces, with 

multiple leading and trailing edge controls, and a fuselage. The mesh generator 

automatically creates the input file for ASP3D by using an extension of the polynomial 

mesh generation procedure developed by Bland.44 Here though, the equations were 

generalized to the physical domain rather than the computational domain, which allows 

the creation of three-dimensional meshes about swept tapered planforms with the ability 

to control the spacing along the leading and trailing edges. This is especially an issue in 

the outboard region of the wing near the tip. For example, Figures 28(a) and 28(b) show 
21

near field views of the planform mesh near the wing tip leading and trailing edges, 

respectively, for the F-5 wing. The figures show that the mesh was generated with a 

uniform spacing about the leading and tailing edges in the streamwise direction, and in 

fact, the spacing is also uniform across the tip in the spanwise direction. The planform 

mesh also emphasizes that all edges of the lifting surface are modeled with grid lines for 

ASP3D application. This is in contrast with meshes generated for CAP-TSD application, 

which cannot have grid lines on leading edges or the edge of the wake (wing tip) because 

of singularities that occur along these lines. Details of the mesh generation procedures 

developed for the ASP3D computer program are presented by Batina.43

To demonstrate the multigrid procedures for a three dimensional case, calculations 

were performed using the 97 x 25 planform mesh (termed Mesh 4) and three subsets of 

the mesh (Meshes 1 through 3) determined by simply deleting every other grid line. This 

is done automatically within the ASP3D code and thus is not something that the user 

needs to do during mesh generation. The four planform meshes are shown in Figures 

29(a) to 29(d) including the total number of points in each direction for each mesh. The 

smallest mesh, for example, is a very coarse 13 x 4 mesh, which has only a few cells on 

the surface of the wing. Similarly, near field views of the corresponding four sectional 

meshes at the root of the F-5 wing are shown in Figures 30(a) through 30(d) with the 

airfoil slit in the center of each plot. The individual figures illustrate the relative mesh 

density near the airfoil surface along the wing root. 

To demonstrate the iterative multigrid capability for time-dependent applications in 

three dimensions, unsteady calculations were performed for the F-5 fighter wing pitching 

about the root midchord with an amplitude of α1 = 0.109° and a reduced frequency based 

on root semichord of k = 0.137. In these calculations the freestream Mach number was

M∞ = 0.899 and the mean angle of attack was α0 = -0.005°. The wing was forced to 

oscillate harmonically in pitch for four cycles of motion as shown in Figure 31(a), to 

ensure periodicity of the flow field, with the resulting lift and moment coefficients shown 

in Figure 31(b). The calculations were performed using only sixteen steps per cycle of 

motion, corresponding to a very large time step of Δt = 1.433, wherein the multigrid 

procedure was used at each time step to reduce the solution residual by two orders of 

magnitude as shown in Figure 31(c). By iterating on the solution at each time step, the 

linearization and factorization errors inherent in the AF2 algorithm are reduced to an 

acceptable level (~10-4) to produce a time accurate solution. The technique is analogous 

to performing subiterations per time step with CAP-TSD to ensure time accuracy, but the 

iterative multigrid procedure in ASP3D is cheaper computationally because the iteration 

is performed on coarser grids within the overall multigrid capability. The corresponding 

number of supersonic points is presented in Figure 31(d). 

The resulting unsteady pressure coefficient distributions, normalized by the 

amplitude of motion and interpolated to the span stations where the NLR experimental 

pressure data40 were measured, are presented in Figure 32. The unsteady pressure 

distributions on the upper surface of the wing are shown in Figure 32(a) and the unsteady 

pressure distributions on the lower surface of the wing are shown in Figure 32(b). The 

ASP3D pressure distributions were computed from the fourth cycle of motion and the 
22

real and imaginary parts of the pressures are the in-phase and out-of-phase components of 

the pressure time histories, respectively. The unsteady pressures on the upper surface of 

the wing agree reasonably well with the experimental data40 except near the shock wave. 

Here the magnitude of the pressure pulse due to the movement of the shock wave is 

slightly over predicted in the ASP3D calculation because of the neglect of viscous 

effects. On the lower surface of the wing (Figure 32(b)) the unsteady pressures are in 

very good agreement with the experimental data, with relatively uniform agreement with 

the data from the wing root to the wing tip. 

Conclusions 

A new computer program has been developed called ASP3D (Advanced Small 

Perturbation – 3D), which solves the small perturbation potential flow equation in an 

advanced form including mass-consistent surface and trailing wake boundary conditions, 

and entropy, vorticity, and viscous effects. The purpose of the program is for unsteady 

aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 

The program exploits the simplicity of stationary Cartesian meshes with the movement or 

deformation of the configuration under consideration incorporated into the solution 

algorithm through a planar surface boundary condition. The ASP3D code is the result of 

a decade of developmental work on improvements to the small perturbation formulation, 

performed while the author was employed as a Senior Research Scientist in the 

Configuration Aerodynamics Branch at the NASA Langley Research Center. The 

ASP3D code is a significant improvement to the state-of-the-art for transonic aeroelastic 

analyses over the CAP-TSD code (Computational Aeroelasticity Program – Transonic 

Small Disturbance), which was developed principally by the author in the mid-1980s. 

The ASP3D code involves a modern cell-centered finite-volume discretization, 

which allows exact planform treatment including leading edges, wing tips, and control 

surface edges. It can treat aircraft configurations involving multiple lifting surfaces with 

leading and trailing edge control surfaces, and a fuselage. The advanced-smallperturbation potential flow equation is solved numerically using either an AF1- or an 

AF2-type approximate factorization algorithm. The AF1 algorithm within the ASP3D 

code is the cell-centered finite-volume version of the AF1 finite-difference algorithm 

developed for CAP-TSD. However, the AF1 algorithm is susceptible to the so-called 

moving shock instability if shock waves move more than one grid point or cell per time 

step. In contrast, the AF2 scheme has a temporal damping term that eliminates the 

moving shock instability by design, and thus, the AF2 algorithm is more robust than the 

AF1 scheme. Thus, with the AF2 approach, convergence to steady state is achieved with 

large time steps or large CFL numbers. And because the moving shock instability does 

not occur, the AF2 algorithm of the ASP3D code is amenable to multigrid 

implementation. Therefore, a multigrid procedure was developed within ASP3D that 

may be used for steady or iterative unsteady applications. The procedure was 

implemented in FAS (full approximation scheme) form, which is applicable to nonlinear 

governing equations. The multigrid procedure may be used for unsteady calculations 

when local iteration is applied. Therefore, the physical time step may be selected for the 

accuracy of the problem under consideration independent of numerical stability issues, 
23

and a large computational step size or CFL number may be selected for rapid 

convergence of the iterative procedure. 

Unsteady aerodynamic and aeroelastic applications of ASP3D were presented to 

assess the time dependent capability and demonstrate various features of the code. The 

cases considered included: (1) the NACA 0012 airfoil undergoing a forced harmonic 

pitching motion about the quarter chord at transonic conditions, (2) a thickening-thinning 

parabolic arc airfoil at a transonic Mach number, (3) an investigation of wave 

propagation characteristics with emphasis on high wave number applications, (4) 

aeroelastic transient calculations for the NACA 0012 airfoil at values of the dynamic 

pressure below, near, and above the flutter value, and (5) the F-5 fighter wing undergoing 

a rigid pitching motion about the wing root midchord axis. The results compared well 

with alternative methods and experimental data, and thus demonstrated the efficiency and 

utility of the ASP3D code for various unsteady aerodynamic and aeroelastic applications. 

Acknowledgments 

The NASA Langley Research Center supported the work contained herein while the 

author was employed as a Senior Research Scientist in the Configuration Aerodynamics 

Branch. The author thanks Dr. David M. Schuster, Head, Structures and Thermal 

Systems Branch, formerly Chief Engineer, Research and Technology Directorate, for 

reviewing the original draft and help in getting the manuscript published. 

References 

1. Bennett, R. M.; Bland, S. R.; Batina, J. T.; Gibbons, M. D.; and Mabey, D. G.: 

Calculation of Steady and Unsteady Pressures on Wings at Supersonic Speeds 

with a Transonic Small-Disturbance Code, AIAA Paper No. 87-0851, April 

1987. 

2. Cunningham, H. J.; Batina, J. T.; and Bennett, R. M.: Modern Wing Flutter 

Analysis by Computational Fluid Dynamics Methods, ASME Paper No. 87-

WA/Aero-9, December 13-18, 1987; also Journal of Aircraft, Vol. 25, October 

1988, pp. 962-968. 

3. Bennett, R. M.; Batina, J. T.; and Bland, S. R.: Wing-Flutter Calculations with 

the CAP-TSD Unsteady Transonic Small-Disturbance Program, AIAA Paper No. 

88-2347, April 1988; also Journal of Aircraft, Vol. 26, September 1989, pp. 876-

882. 

4. Gibbons, M. D.: Aeroelastic Calculations using CFD for a Typical Business Jet 

Model, NASA CR-4753, September 1996. 

5. Batina, J. T.; Seidel, D. A.; Bland, S. R.; and Bennett, R. M.: Unsteady Transonic 

Flow Calculations for Realistic Aircraft Configurations, AIAA Paper No. 87-

0850, April 1987; also Journal of Aircraft, Vol. 26, January 1989, pp. 21-28. 
24

6. Batina, J. T.: An Efficient Algorithm for Solution of the Transonic SmallDisturbance Equation, AIAA Paper No. 87-0109, January 1987; also Journal of 

Aircraft, Vol. 25, July 1988, pp. 598-605. 

7. Batina, J. T.: A Finite-Difference Approximate-Factorization Algorithm for 

Solution of the Transonic Small-Disturbance Equation, NASA TP 3129, January 

1992. 

8. Gibbons, M. D.: personal communication, October 20, 1994. 

9. Steinhoff, J.; and Jameson, A.: Multiple Solutions of the Transonic Potential 

Flow Equation, AIAA Paper No. 81-1019, June 1981. 

10. Williams, M. H.; Bland, S. R.; and Edwards, J. W.: Flow Instabilities in 

Transonic Small Disturbance Theory, NASA TM 86251, January 1985. 

11. Batina, J. T.: Unsteady Transonic Algorithm Improvements for Realistic Aircraft 

Applications, AIAA Paper No. 88-0105, January 1988; also Journal of Aircraft, 

Vol. 26, February 1989, pp. 131-139. 

12. Batina, J. T.: Unsteady Transonic Small-Disturbance Theory Including Entropy 

and Vorticity Effects, AIAA Paper No. 88-2278, April 1988; also Journal of 

Aircraft, Vol. 26, June 1989, pp. 531-538. 

13. Batina, J. T.: Advanced Small Perturbation Potential Flow Theory for Unsteady 

Aerodynamic and Aeroelastic Analyses, NASA/TM-2005-213908, 

November 2005. 

14. Whitfield, D. L.: Analytical Description of the Complete Two-Dimensional 

Turbulent Boundary-Layer Velocity Profile, AEDC-TR-77-79, September 1977; 

also AIAA Paper No. 78-1158, July 1978. 

15. Swafford, T. W.: Analytical Approximation of Two-Dimensional Separated 

Turbulent Boundary-Layer Velocity Profiles, AEDC-TR-79-99, October 1980. 

16. Drela, M.: Two-Dimensional Transonic Aerodynamic Design and Analysis 

Using the Euler Equations, GTL Report No. 187, February 1986. 

17. Batina, J.T.: Introduction to the ASP3D Computer Program for Unsteady 

Aerodynamic and Aeroelastic Analyses, NASA/TM-2005-213909, 

December 2005. 

18. Kuethe, A. M.; and Chow, C.-Y.: Foundations of Aerodynamics: Bases of 

Aerodynamic Design, John Wiley and Sons, New York, Third Edition, 1976. 
25

19. Hafez, M.; and Lovell, D.: Entropy and Vorticity Corrections for Transonic 

Flows, AIAA Paper No. 83-1926, July 1983. 

20. Clebsch, A.: Uber die Integration der Hydrodynamischen Gleichungen, J. Reine 

Angew. Math., Vol. 57, 1959, pp. 1-10. 

21. Grossman, B.: The Computation of Inviscid Rotational Gasdynamic Flows Using 

an Alternate Velocity Decomposition, AIAA Paper No. 83-1900, July 1983. 

22. Dang, T. Q.; and Chen, L. T.: An Euler Correction Method for Two- and ThreeDimensional Transonic Flows, AIAA Paper No. 87-0522, January 1987. 

23. Fuglsang, D. F.; and Williams, M. H.: Non-Isentropic Unsteady Transonic Small 

Disturbance Theory, AIAA Paper 85-0600, April 1985. 

24. Howlett, J. T.: Efficient Self-Consistent Viscous-Inviscid Solutions for Unsteady 

Transonic Flow, AIAA Paper No. 85-0482, January 1985. 

25. Howlett, J. T.: Calculation of Unsteady Transonic Flows With Mild Separation 

by Viscous-Inviscid Interaction, NASA TP 3197, June 1992. 

26. Edwards, J. W.: Transonic Shock Oscillations Calculated with a New Interactive 

Boundary Layer Coupling Method, AIAA Paper No. 93-0777, January 1993. 

27. Godunov, S.: A Finite Difference Method for the Numerical Computation of 

Discontinuous Solutions of the Equations of Fluid Dynamics, Mathematics of the 

USSR, Sbornik, Vol. 47, 1959. 

28. Engquist, B. E.; and Osher, S. J.: Stable and Entropy Satisfying Approximations 

for Transonic Flow Calculations, Mathematics of Computation, Vol. 34, No. 149, 

January 1980, pp. 45-75. 

29. Murman, E. M.: Analysis of Embedded Shock Waves Calculated by Relaxation 

Methods, Proceeds of the AIAA Computational Fluid Dynamics Conference, 

July 1973, pp. 27-40. 

30. Ballhaus, W. F.; Jameson, A.; and Albert, J.: Implicit Approximate-Factorization 

Schemes for Steady Transonic Flow Problems, AIAA Journal, Vol. 16, June 

1978, pp. 573-579. 

31. Jameson, A.: Acceleration of Transonic Potential Flow Calculations on Arbitrary 

Meshes by the Multiple Grid Method, AIAA Paper No. 79-1458, July 1979. 

32. Brandt, A.: Multi-Level Adaptive Solutions to Boundary-Value Problems, 

Mathematics of Computation, Vol. 31, April 1977, pp. 333-390. 
26

33. Anderson, W. K.: Implicit Multigrid Algorithms for the Three-Dimensional Flux 

Split Euler Equations, Doctoral Thesis, Mississippi State University, Mississippi, 

August 1986. 

34. Landon, R. H.: “NACA 0012. Oscillating and Transient Pitching,” Data Set 3 in 

AGARD-R-702, Compendium of Unsteady Aerodynamic Measurements, August 

1982. 

35. Goorjian, P. M.: Implicit Computations of Unsteady Transonic Flow Governed 

by the Full Potential Equation in Conservation Form, AIAA Paper No. 80-0150, 

January 1980. 

36. Tijdeman, H.; and Zwaan, R. J.: On the Prediction of Aerodynamic Loads on 

Oscillating Wings in Transonic Flow, National Aerospace Laboratory, the 

Netherlands, NLR MP 73026 U, April 1973. 

37. Rowe, W. S.; and Cunningham, H. J.: On the Convergence of Unsteady 

Generalized Aerodynamic Forces, Journal of Aircraft, Vol. 21, June 1984, 

pp.420-427. 

38. Robinson, B. A.; Batina, J. T.; and Yang, H. T. Y.: Aeroelastic Analysis of 

Wings Using the Euler Equations with a Deforming Mesh, Journal of Aircraft, 

Vol. 28, November 1991, pp. 778-788. 

39. Isogai, K.; and Suetsuga, K.: Numerical Simulation of Transonic Flutter of a 

Supercritical Wing, National Aerospace Laboratory, Japan, Report TR-276T, 

August 1982. 

40. Tijdeman, H.; Van Nunen, J. W. G.; Kraan, A. N.; Persoon, A. J.; Poestkoke, R.; 

Roos, R.; Schippers, P.; and Siebert, C. M.: Transonic Wind Tunnel Tests on an 

Oscillating Wing with External Stores, AFFDL-TR-78-194, December 1978. 

41. Bland, S. R.: Development of Low-Frequency Kernel-Function Aerodynamics 

for Comparison with Time-Dependent Finite-Difference Methods, NASA TM 

83283, May 1982. 

42. Anderson, W. K.; Thomas, J. L.; and van Leer, B.: Comparison of Finite-Volume 

Flux-Vector Splittings for the Euler Equations, AIAA Paper No. 85-0122, 

January 1985. 

43. Batina, J. T.: Automated Mesh Generation for the ASP3D Computer Program, to 

be published. 

44. Bland, S. R.: Interactive Grid Generation Program for CAP-TSD, NASA TM 

102705, October 1990. 
27

(a) Total mesh. 

(b) Near field view. 

Figure 1 – Cartesian finite-volume mesh (193 x 129) with twenty-chord extent for 

airfoil unsteady aerodynamic and aeroelastic applications with the ASP3D 

computer program. 
28

(a) Forced harmonic pitching for four cycles of motion. 

(b) Lift and moment coefficient responses. 

Figure 2 – Iterative multigrid calculation using 32 steps per cycle of motion for the 

NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 

M∞ = 0.755, α0 = 0.016°, α1 = 2.51°, and k = 0.0814.
29

(c) Last cycle of the iterative multigrid residual history. 

(d) Last cycle of the number of supersonic points.

Figure 2 – Concluded.

iterations 

iterations 

32 steps/cycle 
30

(a) Forced harmonic pitching for four cycles of motion. 

(b) Lift and moment coefficient responses. 

Figure 3 – Iterative multigrid calculation using 16 steps per cycle of motion for the 

NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 

M∞ = 0.755, α0 = 0.016°, α1 = 2.51°, and k = 0.0814.
31

(c) Last cycle of the iterative multigrid residual history. 

(d) Last cycle of the number of supersonic points.

Figure 3 – Concluded.

iterations 

iterations

16 steps/cycle 
32

(a) Forced harmonic pitching for four cycles of motion. 

(b) Lift and moment coefficient responses. 

Figure 4 – Iterative multigrid calculation using 8 steps per cycle of motion for the 

NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 

M∞ = 0.755, α0 = 0.016°, α1 = 2.51°, and k = 0.0814.
33

(c) Last cycle of the iterative multigrid residual history. 

(d) Last cycle of the number of supersonic points.

Figure 4 – Concluded. 

iterations 

iterations 

8 steps/cycle 
34

(a) Forced harmonic pitching for four cycles of motion.

(b) Lift and moment coefficient responses. 

Figure 5 – Iterative multigrid calculation using 4 steps per cycle of motion for the 

NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 

M∞ = 0.755, α0 = 0.016°, α1 = 2.51°, and k = 0.0814.
35

(c) Last cycle of the iterative multigrid residual history. 

(d) Last cycle of the number of supersonic points.

Figure 5 – Concluded. 

iterations

iterations 

4 steps/cycle 
36

(a) Lift coefficient. 

(b) Moment coefficient. 

Figure 6 – Effects of step size on the lift and moment coefficient responses in the last 

(fourth) cycle of motion for the NACA 0012 airfoil oscillating harmonically in pitch 

about the quarter chord at M∞ = 0.755, α0 = 0.016°, α1 = 2.51°, and k = 0.0814.
37

Figure 7 – Definition of midchord thickness as a function of time for the thickeningthinning parabolic arc airfoil.

δ(t) = 0.110 −15 t

15

⎛

⎝

⎜ ⎞

⎠

⎟ + 6 t

15

⎛

⎝

⎜ ⎞

⎠

⎟

⎡ 2

⎣

⎢

⎢

⎤

⎦

⎥

⎥

t

15

⎛

⎝

⎜ ⎞

⎠

⎟

3

for 0 ≤ t ≤15

δ(t) = 0.110 −15 30 − t

15

⎛

⎝

⎜ ⎞

⎠

⎟ + 6 30 − t

15

⎛

⎝

⎜ ⎞

⎠

⎟

⎡ 2

⎣

⎢

⎢

⎤

⎦

⎥

⎥

30 − t

15

⎛

⎝

⎜ ⎞

⎠

⎟

3

for 15 ≤ t ≤ 30

δ(t) = 0 for t ≥ 30

δ(t) 
38

(a) Residual history. 

(b) Number of supersonic points. 

Figure 8 – Iterative multigrid calculation for the thickening-thinning parabolic arc 

airfoil problem at M∞ = 0.85 (32 total time steps). 

iterative 

multigrid 
39

(a) Residual history. 

(b) Number of supersonic points. 

Figure 9 – Iterative unigrid calculation for the thickening-thinning parabolic arc 

airfoil problem at M∞ = 0.85 (32 total time steps).

iterative 

unigrid 
40

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 10 – Instantaneous results at t = 8.5 from the iterative multigrid 

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85.

t = 8.5 
41

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 11 – Instantaneous results at t = 11.5 from the iterative multigrid 

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85.

t = 11.5 
42

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 12 – Instantaneous results at t = 18.25 from the iterative multigrid

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85. 

t = 18.25 
43

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 13 – Instantaneous results at t = 26.875 from the iterative multigrid 

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85.

t = 26.875 
44

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 14 – Instantaneous results at t = 29.0 from the iterative multigrid 

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85.

t = 29.0 
45

(a) Pressure coefficient distribution. 

(b) Near field pressure coefficient contour lines. 

Figure 15 – Instantaneous results at t = 32.0 from the iterative multigrid 

calculation of the thickening-thinning parabolic arc airfoil at M∞ = 0.85.

t = 32.0 
46

(a) ASP3D pressure coefficient distributions at t = 8.5, 11.5, and 18.25. 

(b) Corresponding FP and TSD pressure coefficients reported by Goorjian.35

Figure 16 – Pressure coefficient comparisons for the thickening-thinning parabolic 

arc airfoil at M∞ = 0.85. 
47

(c) ASP3D pressure coefficient distributions at t = 26.875, 29.0, and 32.0. 

(d) Corresponding FP and TSD pressure coefficients reported by Goorjian.35

Figure 16 – Concluded.
48

(a) Subsonic freestream Mach number of M∞ = 0.8. 

(b) Transonic freestream Mach number of M∞ = 0.875. 

Figure 17 – NLR diagrams36 of upstream propagation of disturbance wave fronts 

generated by a source at the control surface hinge line. 
49

(a) Subsonic freestream Mach number of M∞ = 0.8.

(b) Transonic freestream Mach number of M∞ = 0.875.

Figure 18 – NLR shadowgraph (spark exposure) photographs36 showing upstream 

wave propagation. 
50

(a) Upstream propagation. 

(b) Downstream propagation. 

Figure 19 – Linear theory wave propagation characteristics as a function of Mach 

number for a range of values of reduced frequency k.
51

(a) Upstream propagation. 

(b) Downstream propagation. 

Figure 20 – Linear theory wave propagation characteristics as a function of Mach 

number for various values of wavelength λ.
52

(a) Real part. 

(b) Imaginary part. 

Figure 21 - Effects of Mach number on linear theory lifting pressure coefficient 

distributions for a flat plate airfoil oscillating about the leading edge at k = 0.9. 
53

(a) First-order-accurate φxt term.

(b) Second-order-accurate φxt term. 

Figure 22 – Effects of order of accuracy on ASP3D lifting pressure coefficient 

distributions for a flat plate airfoil oscillating about the l. e. at k = 0.9 and M∞ = 0.9. 
54

(a) ASP3D pressure distribution including entropy and vorticity effects. 

(b) Euler (CFL3D) pressure distribution computed by Robinson, et al.38

Figure 23 – Steady pressure coefficient distributions for the NACA 0012 airfoil at 

M∞ = 0.8 and α0 = 0°. 

ASP3D 
55

Figure 24 – Generalized displacements corresponding to stable (converging) 

aeroelastic motion for the NACA 0012 airfoil at M∞ = 0.8, α0 = 0°, and Q = 0.2 for 

the structural parameter values of Isogai Case “A”. 

Q = 0.2 
56

Figure 25 – Generalized displacements corresponding to near-neutrally stable 

(flutter) aeroelastic motion for the NACA 0012 airfoil at M∞ = 0.8, α0 = 0°, and 

Q = 0.5 for the structural parameter values of Isogai Case “A”.

Q = 0.5 
57

Figure 26 – Generalized displacements corresponding to unstable (diverging) 

aeroelastic motion for the NACA 0012 airfoil at M∞ = 0.8, α0 = 0°, and Q = 0.8 for 

the structural parameter values of Isogai Case “A”. 

Q = 0.8
58

(a) Near field view of 97 x 25 planform mesh. 

(b) Near field view of 97 x 65 root sectional mesh. 

Figure 27 – Finite volume meshes for the F-5 fighter wing.
59

(a) Wing tip leading edge region.

(b) Wing tip trailing edge region. 

Figure 28 – Near field views of the planform mesh of the F-5 fighter wing.

F-5 

wing

wing 

 tip

F-5 

wing

wing 

 tip
60

(a) Mesh 4 (97 x 25). (b) Mesh 3 (49 x 13). 

(c) Mesh 2 (25 x 7). (d) Mesh 1 (13 x 4). 

Figure 29 – Near field view of finite volume planform meshes for the F-5 fighter 

wing for ASP3D multigrid calculations.
61

(a) Mesh 4 (97 x 65). (b) Mesh 3 (49 x 33). 

(c) Mesh 2 (25 x 17). (d) Mesh 1 (13 x 9). 

Figure 30 – Near field view of finite volume root sectional meshes for the F-5 fighter 

wing for ASP3D multigrid calculations. 
62

(a) Forced harmonic pitching for four cycles of motion. 

(b) Lift and moment coefficient responses. 

Figure 31 – Iterative multigrid calculation using 16 steps per cycle of motion for the 

F-5 fighter wing pitching about the root midchord at M∞ = 0.899, 

α0 = -0.005°, α1 = 0.109°, and k = 0.137.
63

(c) Iterative multigrid residual history.

(d) Number of supersonic points.

Figure 31 – Concluded. 

iterations 

iterations 

16 steps/cycle 
64

(a) Upper surface. 

(b) Lower surface. 

Figure 32 - ASP3D inviscid unsteady pressure coefficient comparisons with NLR 

experimental data40 for the F-5 fighter wing pitching at M∞ = 0.899, 

α0 = -0.005°, α1 = 0.109°, and k = 0.137.

 F-5 

wing

 F-5 

wing 

M∞ = 0.899

M∞ = 0.899

 ASP3D inviscid solution 

 NLR data - real part 

 NLR data - imag. part

 ASP3D inviscid solution 

 NLR data - real part 

 NLR data - imag. part
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Application of the ASP3D Computer Program to Unsteady Aerodynamic

and Aeroelastic Analyses

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Batina, John T.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-19234

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

An electronic version can be found at http://ntrs.nasa.gov

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 02

Availability: NASA CASI (301) 621-0390

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

A new computer program has been developed called ASP3D (Advanced Small Perturbation - 3D), which solves the small

perturbation potential flow equation in an advanced form including mass-consistent surface and trailing wake boundary

conditions, and entropy, vorticity, and viscous effects. The purpose of the program is for unsteady aerodynamic and

aeroelastic analyses, especially in the nonlinear transonic flight regime. The program exploits the simplicity of stationary

Cartesian meshes with the movement or deformation of the configuration under consideration incorporated into the solution

algorithm through a planar surface boundary condition. The paper presents unsteady aerodynamic and aeroelastic applications

of ASP3D to assess the time dependent capability and demonstrate various features of the code.

15. SUBJECT TERMS

Aeroelasticity; Computational Fluid Dynamics; Unsteady Aerodynamics; Computer Programs

18. NUMBER

 OF 

 PAGES

69

19b. TELEPHONE NUMBER (Include area code)

(301) 621-0390

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

843515.02.02.29.07

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2006-214277

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 02 - 2006