Sample Cover
AIAA 2003-0774
Design-Oriented Aerodynamic Analysis
for Supersonic Laminar Flow Wings
I. Kroo, P. Sturdza
Stanford University
Stanford, CA
41st Aerospace Sciences Meeting & Exhibit
6 - 9 January 2003
Reno, Nevada
For permission to copy or to republish, contact the American Institute of Aeronautics and Astronautics,
1801 Alexander Bell Drive, Suite 500, Reston, VA, 20191-4344.
AIAA-2003-0774
Design-Oriented Aerodynamic Analysis for Supersonic Laminar Flow Wings
Ilan Kroo∗
Peter Sturdza†
Department of Aeronautics and Astronautics
Stanford University, Stanford, CA
Abstract
Introduction
The computation of boundary layer properties and
laminar-to-turbulent transition location is a very
complex problem, generally not undertaken in the
context of aircraft design. Yet if an aircraft is to
exploit the advantages of laminar flow and proper
trade-offs between inviscid drag, structural weight,
and skin friction are to be made, this is just what
must be done. This paper summarizes a designoriented method for the aerodynamic analysis of supersonic wings including approximate means for estimating transition and total drag. The boundary
layer analyses employed here are meant to be inexpensive but sufficiently accurate to provide some
guidance for advanced design studies or to be incorporated in a multidisciplinary optimization.
The boundary layer analysis consists of a quasi3D finite difference method, based on inputs from
a 3D Euler analysis. A parametric fit of integral
boundary layer properties to the amplification rate
predicted by linear stability analysis, leads to a simplified method for predicting TS transition and the
onset of instability of stationary crossflow modes.
This simplifed boundary layer analysis and inviscid CFD were combined to evaluate the drag of a
simple wing planform, defined by approximately ten
geometric parameters. These parameters were used
as design variables in a numerical optimization procedure to minimize total wing drag, illustrating the
trade-offs between skin friction and inviscid drag.
Three dimensional Computational Fluid Dynamics
(CFD) has become a useful and popular tool for the
analysis and design of aircraft. However, for multidisciplinary optimization (MDO) at the conceptual
design stage, for cases in which laminar to turbulent
transition play an important role in configuration
drag, it is still far too expensive.
Theory and limited experimental evidence suggest
the possibility of maintaining extensive natural laminar flow over wings with moderate sweep at supersonic speeds. The wing planform, twist, and airfoil
profiles must carefully designed, taking into account
the effects on both aerodynamics and other disciplines. The tradeoffs between different disciplines
are a key element to the viability of the concept.
For example, wings that support the maximum extent of natural laminar flow may have little sweep
and thus require very small thickness to reduce wave
drag. This adversely affects the structural design,
fuel volume, and subsonic performance. A detailed
and complete multidisciplinary design study is necessary to determine whether these tradeoffs result
in a more efficient configuration than a conventional
design. (For a discussion of the MDO problem, see
Manning and Kroo.1 )
Since 3D CFD is incapable of the quick execution
times necessary for this type of design problem, a
less expensive, yet reasonably accurate method of
boundary layer computation within an optimization
framework is sought. Such a method must provide a calculation of skin friction, including sufficient modeling to enable boundary layer transition
prediction, while demanding little computational resources, since the analysis will likely be called thousands of times in the course of a nonlinear search
procedure.
In the following sections, we describe this approach to the design of supersonic wings with lam-
∗ Professor,
† Graduate
Fellow AIAA
Student
c
Copyright 2003
by the authors. Published by the American
Institute of Aeronautics and Astronautics, Inc., with permission.
1
American Institute of Aeronautics and Astronautics
2
inar flow, including boundary layer calculation and
approximate transition prediction. Comparisons
with more refined calculations and application to a
simple wing design problem illustrate the use of the
method.
Direct Optimization Approach
Previous laminar design studies have often relied on
the specification of a target pressure distribution,
followed by nonlinear inverse design to determine a
geometry compatible with laminar flow. The desired Cp may be defined based on the experience
of aerodynamics experts and designers who understand some of the basic trade-offs between inviscid
and viscous drag as well as the need for maintaining certain structural properties. However, recent
experiences with the design of low boom, small supersonic aircraft suggest that this approach quickly
becomes intractable when complex trades between
supersonic and subsonic performance, boom signature shape, and wing structural weight are required.
Since wing skin friction may account for less than
20% of the drag of a business jet design at Mach 1.6,
the implications of the target Cp on properties other
than laminar flow extent may be paramount. Furthermore, if only the outer wing panel is amenable
to extensive laminar flow and/or it is only practical to exploit laminar flow on the upper (or lower)
surface, the importance of laminar flow is further reduced and it is not clear that the best design should
even try to obtain low skin friction at the expense
of other characteristics.
A quantifiable evaluation of this trade-off is possible using a direct optimization approach that involves varying the wing geometric properties and
computing the effect on overall aircraft performance
and constraints. This multidisciplinary optimization approach addresses the difficulties in identifying suitable target pressures and can provide important insight into basic configuration choices. It
has not, however, been used for advanced design of
laminar aircraft for several reasons. The computation of boundary layer properties and transition is
extremely computation intensive and often not reliable. Results may also be rather noisy, causing
potential difficulties with conventional search methods. The goal of the present research was to develop
methods, which although necessarily approximate,
would provide more direct guidance in the design of
supersonic vehicles employing laminar flow for improved performance.
Sweep/Taper Boundary Layer
Analysis
The approach to estimating boundary layer properties is based on an approximate sweep/taper method
described by the authors in.2 Some of the discussion
here is excerpted from that paper, while a more detailed description of the computations are provided
in Sturdza’s thesis.3 The sweep/taper method was
originally suggested by Bradshaw as an extension to
the ideas used for calculating boundary layers on infinite yawed wings. The sweep/taper simplification
allows the computations to be two dimensional even
though a three dimensional flow field is solved. It
applies to semi-infinite straight tapered, swept wings
with straight isobars parallel to the spanwise generator lines. Bradshaw et al.,4, 5 apply this simplification to their turbulent boundary layer code and
report that it runs at approximately two thirds the
speed of their compressible 2D code. It does not
compute laminar boundary layers due to the special
numerical procedure of the Bradshaw-Ferriss-Atwell
finite difference method. Kaups and Cebeci6 developed a laminar finite difference boundary layer program based on similar assumptions, however, it can
only calculate the boundary layer up to the chordwise station where the crossflow changes sign (unacceptable for the current application). In the mideighties, Ashill and Smith7 developed an extension
to Smith’s general 3D integral method that takes advantage of the Bradshaw sweep/taper theory. This
approach, however, introduces additional assumptions that also make it unsuitable for this research.
The main difficulty is, again, an inability to handle configurations where the crossflow changes sign.
The present method is capable of laminar and turbulent calculations with a slightly improved version of
Bradshaw’s formulation of the sweep/taper theory.
The basic sweep/taper theory can be developed
as follows. In the case of the infinite (untapered)
yawed wing for incompressible laminar flow, it is
well known that the boundary layer equations become two dimensional in the plane normal to the
spanwise coordinate, with a third, decoupled equation for the spanwise velocity (Prandtl’s indepen-
3
∂u ∂v
+
=0
∂x ∂y
∂
∂
∂p
∂
∂u
x-mom: ρ u
+v
u=−
+
µ
∂x
∂y
∂x ∂y
∂y
∂
∂
∂
∂w
z-mom: ρ u
+v
w=
µ
,
∂x
∂y
∂y
∂y
U
8
dence principle):
cont:
where (u, v, w) and (x, y, z) are the cartesian velocity
components and coordinates in the streamwise, normal, and spanwise directions, respectively. In the
compressible or turbulent case, the equations become coupled, but since the derivatives in the spanwise direction are zero, the equations are still computationaly two dimensional (there are only two independent variables). After adding taper to the infinite yawed wing, as discovered by Bradshaw et al.,4
the boundary layer equations can be transformed
in order to remain computationaly two dimensional
(with the additional assumption that the surface isobars are parallel to the spanwise generators of the
wing). The 3D compressible boundary layer equations are written in cylindrical coordinates along an
arc from the leading edge to the trailing edge that
is perpendicular to both edges and to the isobars in
between (Figure 1):
cont:
x-mom:
z-mom:
temp:
∂ρu ∂ρv ∂ρw ρw
+
+
−
=0
∂x
∂y
∂z
r
o
∂
∂
∂
ρuw
ρ u
+v
+w
u−
∂x
∂y
∂z
ro
∂p
∂
∂u
=−
+
µ
∂x ∂y
∂y
∂
∂
∂
ρu2
ρ u
+v
+w
w+
∂x
∂y
∂z
ro
∂p
∂
∂w
=−
+
µ
∂z
∂y
∂y
∂
∂
∂
ρcp u
+v
+w
T
∂x
∂y
∂z
∂p
∂p
∂
∂T
=u
+w
+
κ
∂x
∂z
∂y
∂y
" 2 #
2
∂u
∂w
+µ
+
.
∂y
∂y
The origin of the coordinate system is where the
wing would come to a point if it were semi-infinite
in span. In this coordinate system, there is no pressure gradient in the spanwise, z, direction, and other
derivatives in the z direction can be written in terms
of derivatives in the direction normal to the wing
z
r0
Λ
x
θ
Figure 1: Cylindrical coordinate system used in
the sweep/taper boundary layer calculation (view of
wing planform from above).
surface which results in a set of computationally 2D
equations. This comes from noting that along a generator line (constant x), the boundary layer thickness is closely proportional to r0 in turbulent flow,4
√
and r0 in laminar flow; essentially a conical flow
assumption where the boundary layer can be viewed
as self similar in the spanwise direction. Therefore
the spanwise derivative, ∂f /∂z, of any flow quantity, f , (other than pressure) can be re-written in
terms of a vertical derivative of that same flow quantity, (y/r0 )∂f /∂y for turbulent flow or (y/2r0 )∂f /∂y
for laminar flow. These operations on the boundary layer equations result in numerically two dimensional equations that can be solved on a 2D x − y
grid along an arc of constant radius from the origin
of the coordinate system.
The Bradshaw sweep/taper assumption limits,
strictly speaking, the isobars to be parallel with the
spanwise generators and thus constrains the lift distribution to a triangle. For the semi-infinite tapered
wing, calculations along a single arc will provide the
3D flow solution on the entire wing. However, by
running calculations at several spanwise stations and
including some corrections to the basic theory described below, this approach can be used as a quasi3D calculation on a finite wing, taking into account
the local effects of sweep and taper as long as the
isobars are reasonably straight and parallel with the
generators. Fortunately, it is not uncommon for the
isobars to be somewhat parallel with the generator lines in both subsonic and supersonic wings of
medium to high aspect ratio.
The present code is an extension of the Bradshaw
sweep/taper theory. Bradshaw’s method, as has
been implemented previously, only uses freestream
conditions and surface pressures along the computational arc as inputs. The external (inviscid) veloc-
4
ity components are calculated from those inputs using the sweep/taper assumptions. It has been found
that the method can be improved by using more
information from the inviscid solution. The entire
inviscid flow field is available as an external boundary condition to the sweep/taper code. For example,
the pressure gradient in the generator direction, ∂p
∂z ,
which is neglected in the original method, is calculated by the inviscid flow solver and can be applied
directly in the boundary layer equations. The external velocity components are similarly available, but
the difficulty is that they are not necessarily consistent with the sweep/taper equations. That deviation, however, can be used to develop corrections to
the sweep/taper version of the mass, momentum and
energy equations in the form of source terms. The
present sweep/taper method implements these improvements by solving the following set of equations
for laminar flow:
∂ρu ∂ρˆ
v 3ρw
+
−
= ρQ
∂x
∂y
2ro
∂
∂
ρuw
x-mom:
ρ u
+ vˆ
u−
∂x
∂y
ro
∂
∂u
∂p
+
µ
+ Fx
=−
∂x ∂y
∂y
∂
ρu2
∂
z-mom:
ρ u
+ vˆ
w+
∂x
∂y
ro
∂p
∂
∂w
+
µ
+ Fz
=−
∂z
∂y
∂y
∂
∂
temp: ρcp u
+ vˆ
T
∂x
∂y
∂p
∂
∂T
∂p
+w
+
κ
=u
∂x
∂z
∂y
∂y
" 2 #
2
∂u
∂w
+µ
+
+ B,
∂y
∂y
cont:
where,
vˆ = v + w
y
,
2ro
and,
∂ρw
y ∂ρw
−
∂z
2ro ∂y
∂u
y ∂u
Fx = −ρw
−
∂z
2ro ∂y
∂w
y ∂w
Fz = −ρw
−
∂z
2ro ∂y
∂T
y ∂T
B = −ρcp w
−
.
∂z
2ro ∂y
Q=−
1
ρ
(1)
(2)
(3)
(4)
Mach
alt.
ΛLE
ΛT E
b/2
cr
ct
t/c
(ft.)
(deg.)
(deg.)
(ft.)
(ft.)
(ft.)
(%)
Agrawal
2.0
40,000
32.0
-15.4
18.94
22.5
5.45
5.0
NLF
1.5
45,000
16.7
-29.0
23.15
25.3
5.52
1.75
Flight Test
1.4
10,000
19.5
-27.5
3.0
4.16
1.57
2.5
Table 1: Geometry and flight conditions for the calculations appearing in this paper.
The source terms, equations 1-4, and the spanwise pressure gradient, ∂p
∂z , are the corrections applied to the Bradshaw sweep/taper theory to arrive
at the present formulation. The pressure gradient is
independent of y, and is therefore known throughout the boundary layer, but the source terms are
not. At the edge of the boundary layer, the source
terms can be evaluated using the inviscid flow solution from an Euler or panel code. At the wing
surface, w and y are zero, so the terms vanish there.
In this implementation, they are modeled linearly in
between.∗ These developments significantly improve
sweep/taper results without any discernible computational cost over the original Bradshaw formulation
of the theory.
The sweep/taper equations are discretized and
solved using Keller and Cebeci’s algorithm for the
compressible boundary layer equations.8 In brief,
the code solves the continuity, x- and z-momentum,
and temperature equations in cylindrical coordinates
along the circular arc drawn in Figure 1. The continuity equation is integrated out with stream functions and the remaining equations are transformed
via a compressible Falkner-Skan procedure using the
arc length and the tangential external velocity component in the similarity variables. The discretization is second order, and Newton’s method is used
to march the solution along the arc from leading
edge to trailing edge. Theoretically, the only limit
on wing sweep and taper is that the marching always
be downwind.
Figure 2 presents the crossflow Reynolds number predicted using the sweep/taper method compared to detailed 3D Navier-Stokes computations
(CFL3D)9 on the Agrawal wing planform (Table 1).
Inviscid pressures and velocities used for all the
boundary layer work in this paper were performed
∗ Actually, they are linear in the transformed Falkner-Skan
vertical coordinate.
5
1
1
Present Sweep/Taper
3D Navier-Stokes
Present Sweep/Taper
3D Navier-Stokes
0.8
0.6
0.6
y/δ
Sweep/Taper without corrections
Sweep/Taper with corrections
3D Navier-Stokes
1000
0.8
y/δ
1200
800
0.4
0.4
0.2
0.2
0
-0.025
Rcf
0
0
600
0.2
0.4
0.6
0.8
1
u/ue
400
0.6
y/δ
y/δ
0.6
0.4
0.4
0.2
0.4
0.6
0.8
0
0.8
0.8
0
-0.005
Present Sweep/Taper
3D Navier-Stokes
Present Sweep/Taper
3D Navier-Stokes
0
-0.015
-0.01
w/ue
1
1
200
-0.02
1
X/C
0.2
0.2
0
0
0.2
0.4
0.6
u/ue
Figure 2: Crossflow Reynolds number: sweep/taper
vs. 3D Navier-Stokes on the Agrawal wing at midsemispan.
Laminar
3D Navier-Stokes
Euler + Sweep/Taper
Skin Friction
0.00043
0.00043
Inviscid
0.01923
0.01919
Turbulent
3D Navier-Stokes
Euler + Sweep/Taper
Skin Friction
0.00362
0.00356
Inviscid
0.01931
0.01919
Table 2: Wing drag coefficients: Sweep/taper vs.
Navier-Stokes on the Agrawal wing for 100% laminar
and 100% turbulent computations.
with the FLO87 Euler code developed by Jameson.10
The present sweep/taper calculations are in quite
good agreement with Navier-Stokes, showing the
dramatic improvement over the original Bradshaw
formulation. Not shown are Agrawal’s results11
which cannot go beyond the zero-sweep chord location because he used a modification of the Kaups
and Cebeci code. The velocity profiles calculated
near mid-chord and near the trailing edge with the
present sweep/taper method also compare favorably
with Navier-Stokes results (Figure 3). It has been
found that subtleties in the flow are captured at
least qualitatively by this method. For instance,
variations at different spanwise stations (such as
wing root and tip effects—departures from strict
sweep/taper assumptions) are computed accurately.
This is shown in Figure 4 on an example NLF
planform (Table 1), verifying that the sweep/taper
method can be used as a “quasi-3D” theory by per-
0.8
1
0
-0.006
0
0.006
0.012
w/ue
Figure 3: Streamwise (left) and crossflow (right) velocity profiles at 56% chord (top) and 92% chord
(bottom) calculated with sweep/taper theory (solid
line) vs. 3D Navier-Stokes (symbols) on the Agrawal
wing.
forming calculations at various spanwise locations.
The coefficient of friction in the freestream direction is compared between the sweep/taper method
and 3D Navier-Stokes in Figure 5 at one spanwise
station. Complete wing drag coefficients are shown
in Table 2 for laminar and turbulent computations
on the Agrawal wing. The agreement is very good
especially in the laminar flow case. Since different
turbulence models are involved (Spalart-Almuras in
CFL3D and Cebeci-Smith in the sweep/taper code)
the difference of less than one count in turbulent skin
friction is very reasonable. The comparisons confirm
that the sweep/taper approach can accurately calculate friction drag on thin supersonic wings and that
there is no need for viscous-inviscid iteration at these
flight conditions.
Figure 6 depicts the temperature profile in the
boundary layer calculated with the sweep/taper
method and with Navier-Stokes. The two sets of
Navier-Stokes data are for the same run, but at two
convergence levels (at overall residual levels of 10−6
and 10−8 ). It takes an order of magnitude more
iterations to achieve temperature convergence compared with the ususal CL or CD convergence criteria
— two processor-weeks on a 195 MHz SGI Origin in
this case. This behavior has consistently been observed in 2D Navier-Stokes calculations as well.12
6
400
350
Sweep/Taper Outboard
Sweep/Taper Inboard
3D Navier−Stokes Outboard
3D Navier−Stokes Inboard
10-2
SWPTPRBL
CFL3D
250
200
150
Cf
Crossflow Reynolds Number
300
10-3
100
50
0
0
0.1
0.2
0.3
0.4
0.5
X/C
0.6
0.7
0.8
0.9
1
10-4
0
Figure 4: Crossflow Reynolds number: sweep/taper
(lines) vs. 3D Navier-Stokes (symbols) at two spanwise stations on an NLF wing.
Since temperature strongly influences the stability
of Tollmien-Schlichting waves, full Navier-Stokes solutions must be highly converged to be useful in stability calculations. For this reason, the speed advantage of the present sweep/taper method over 3D
Navier-Stokes is quite large. A week or two of runtime can be shortened to about 12 seconds on the
same computer.
Transition Prediction
Transition to turbulence is one of the most complex
aspects of fluid dynamics, however, reasonable semiempirical methods exist that can capture the general
behavior of the transition front. One of the most
common of such methods is linear stability theory
and the en transition criteria. The en method is the
basis of the current work. As additional 3D supersonic experimental data become available, improved
techniques can be substituted in the general framework developed here.
Transition to turbulence in a boundary layer
is caused by growing disturbances in the laminar
boundary layer. These flow disturbances are initiated mainly by surface roughness or undulations,
mechanical and acoustic vibrations and by shock
waves or turbulence in the free-stream flow. The
disturbances start out as small, periodic perturbations to the mean flow, then grow or decay depending upon the characteristics of the boundary-layer
velocity and temperature profiles. At first, when the
disturbances are small, their behavior is similar to
sinusoidal plane waves. As the unstable modes grow
0.2
0.4
0.6
0.8
1
X/C
Figure 5: Skin friction coefficient: Sweep/taper vs.
Navier-Stokes on the Agrawal wing.
in amplitude, nonlinear interactions become dominant. Following the nonlinear growth, these laminar
disturbances begin to spurt out intermittent spots
of turbulence which spread and eventually merge together resulting in a fully turbulent boundary layer.
A good discussion of the physics of transition can be
found in White.13
These laminar instabilities are classified into several different types including Tollmien-Schlichting
waves (TS), G¨ortler vortices, attachment line instabilities, and crossflow vortices (CF). The focus of
the present work is on the development of approximate methods by which TS and CF transition can
be estimated.
Transition Analysis Overview
A complete analysis of transition can really only
be done with direct numerical simulation (DNS),
however, as in the computation of turbulence, this
method is limited to small Reynolds numbers due to
the enormous computational requirements. Even so,
knowledge of the factors that originally excite these
instabilities is incomplete — the initial or boundary
conditions for the DNS. In the case of an airplane
wing, it is difficult to specify a complete picture of
the conditions: free-stream turbulence, acoustic energy both in the flow and in the structure, microscopic details of the surface shape and surface finish
and contaminants. Receptivity, the response of laminar modes to various disturbances, is an active area
7
1
0.8
Present Sweep/Taper
3D Navier-Stokes (Resid: 10-8)
3D Navier-Stokes (Resid: 10-6)
y/δ
0.6
0.4
0.2
0
200 220 240 260 280 300 320 340 360
T (Kelvin)
Figure 6: Temperature profile showing slow convergence of N-S calculation (the more converged case
required an order of magnitude more iterations).
of research. Although rules of thumb and empirical results from receptivity analysis and experiments
can be very useful, detailed receptivity analysis is
not reasonable at the conceptual and preliminary
design stages.
The solution of stability equations represents a
step down from the brute force DNS approach.
These equations are a form of the Navier-Stokes
equations after subtracting out the steady mean flow
terms, leaving just the description of the unsteady
disturbances. The significant advantage is that with
relatively few assumptions, these equations can be
manipulated to allow very efficient numerical solutions. Since one of the approximations results in
partial differential equations with parabolic character, the equations are known as the parabolized
stability equations (PSE). The simplifying assumptions mainly involve a description of the perturbations to the mean flow as a periodic phenomena
characterized by a slowly varying amplitude function
and rapidly varying wave behavior. The basic idea
is that the oscillatory portions are transformed by
Fourier series, leaving just the amplitude function to
solve. Since the amplitude, or shape function, varies
slowly in the streamwise direction, much less computational resolution is required to capture an accurate
solution compared with DNS solutions. Additionally, a careful choice for the decomposition of the
solution into the wave and shape functions allows
some terms to be neglected, resulting in parabolic
equations that can be solved with very efficient space
marching schemes. Further details can be found in
papers by Malik, Chang and Li.14, 15 The accuracy
of the solution is excellent because laminar instabilities conform very nicely to the simplifying assumptions in the PSE, especially the nonlinear form of
the PSE.16
Although the PSE provide a practical method
with today’s computers to calculate laminar instabilities, there are drawbacks. For transition caused by
streamwise Tollmien-Schlichting waves it is widely
accepted that PSE is not necessary since Linear
Stability Theory (LST, described below) provides
good transition prediction capabilities. For the other
dominant mode in swept wings in atmospheric flight
– stationary crossflow vortices – it is thought that
LST is adequate if the transition amplitude is empirically found for similar geometries.17 Additionally,
it is not simple to predict transition for the stationary crossflow mode using PSE because the crossflow
vortices themselves do not cause transition directly,
but rather saturate and allow a secondary instability
to grow which then triggers transition. This means
that yet another stability code is used to calculate
the growth of this secondary instability and an nfactor correlation is then applied. Another drawback is that the nonlinear PSE require as inputs initial conditions describing the birth of the laminar
instabilities which are not generally available. Frequently, when detailed experimental flow fields are
not known, the PSE solutions are started from mode
shapes derived using LST.
LST, on the other hand, is very well suited for
this work. It can be used as a self-contained, though
approximate, transition prediction method without
requiring extensive knowledge of initial conditions.
The main drawback is that the numerical procedures
are poorly-suited to optimization. This is because
LST solvers need to be run interactively, with an expert user monitoring results, discarding non-physical
solutions, modifying inputs, and helping the solver
along in some cases. Otherwise, LST solvers are reasonably fast, can accommodate the high Reynolds
number flows in this study and do not require special starting procedures. It is true that LST does
not provide as detailed a description of the physical flow disturbances as PSE, however that is not
always an advantage for transition prediction at the
conceptual design stage. Greater resolution of the
flow requires a deeper understanding of the physics
involved in the process of receptivity and transition
itself. The simpler theory can be more robust.
8
LST can be summarized as follows. The NavierStokes equations are the starting point. They are
separated into mean flow terms and fluctuation
terms which are then simplified by neglecting nonlinear effects. Then the mean flow equations are
subtracted out, leaving linearized differential equations governing the small disturbances. The parallel
flow assumption is applied, which is basically a statement that the boundary-layer growth is small over
a wavelength of the disturbance. With the additional assumption of sinusoidal solutions (via complex exponentials), the equations simplify into an
eigenvalue problem at each chordwise station. This
means that the solution procedure involves marching from leading to trailing edge and solving for the
local disturbance properties. The inputs are the disturbance frequency and wave angle, and the outputs
are the amplification rate and mode shapes for those
inputs. Mathematical detail can be found in a paper
by L. M. Mack.18
The difficulty is in the numerical procedure that
finds unstable modes and tracks them as they travel
downstream. This quality makes linear stability theory (LST) codes necessarily interactive, with the
user providing reasonable frequencies and checking
for non-physical results, or modes that are lost from
station to station. For this reason, even a very powerful computer will still not enable the use of LST
codes in numerical design optimization.
Linear stability theory, and for that matter PSE
calculations as well, don’t predict transition, but
rather just predict the growth of periodic laminar
instabilities. The most popular transition criterion
based on LST calculations originated in the mid
1950s. Smith and Gamberoni, and, independently,
van Ingen compared linear stability calculations
with experimental results of various boundary layers
where the transition location was known.19, 20 They
found that the transition location in many experiments coincided with the location where linear stability theory predicts the amplification of TollmienSchlichting waves to reach about 8100. Since the
commonly used amplification n factor is the natural logarithm of the amplification ratio, this became
known as the e9 method.
This method works very well for two-dimensional
boundary layers that experience transition caused
by growth of Tollmien-Schlichting instabilities. It
is also assumed that the free-stream turbulence and
acoustic disturbances are not too great since a long,
slow, linear growth of the instabilities is required for
LST to be accurate. If the initiating disturbances
are large, then nonlinear growth can dominate the
transitional region and violate the assumptions of
LST. Low turbulence tunnels and atmospheric flight
conditions generally meet these requirements. The
n factor can be adjusted to other than 9 to account
for varying receptivity effects, such as free-stream
turbulence or surface roughness. L. M. Mack suggested this approach by proposing a correlation for
n as a function of free-stream turbulence. He argued
that the free-stream turbulence influences the initial
amplitude of the laminar instabilities.21 As a rough
guide, transition in wind tunnels generally observed
at n factors of 7 to 9 whereas in atmospheric flight,
it is more common to see transition at n factors of
10 to 14. For this reason, this criterion is now more
commonly known as the en method where n is found
empirically for the conditions of interest.
Present Transition Prediction Methodology
The present development assumes that 3D transition to turbulence on swept wings can be split into
streamwise and crossflow criteria that are calculated
separately and subsequently recombined. Compressible linear stability theory is heavily used since it
successfully predicts TS dominated transition, and,
it is believed, stationary crossflow transition on
swept wings. Since the goal is to provide a transition analysis that can be used for aircraft design optimization, linear stability results are predetermined
and stored as parametric fits. These fits are evaluated during the optimization to calculate instability
amplification. The en criterion is used to predict
transition for both the streamwise TS calculation
and the crossflow modes (but with different values
of n).
Streamwise Instability Calculation
Figure 7: Boundary-layer instabilities for similarity
solutions.
The streamwise methodology is borrowed from the
work of Drela,22 Gleyzes23 and Stock and Dagen-
9
hart.? Most of these previous works are only applicable to incompressible flows, with Drela extending
the method to the transonic regime.
8
7
6
Present Fit
2500 Hz
2750 Hz
3000 Hz
3250 Hz
3500 Hz
4000 Hz
4500 Hz
5000 Hz
5500 Hz
6000 Hz
7000 Hz
8000 Hz
9000 Hz
10000 Hz
12000 Hz
15000 Hz
17500 Hz
20000 Hz
25000 Hz
30000 Hz
5
15
N
M=1.5 t/c=1%
4
Amplification Ratio, n
3
2
M=1.5 t/c=2%
10
1
M=2.4 t/c=1%
M=1.5 t/c=3%
0
0.5
1
1.5
2
2.5
x (m)
5
M=2.4 t/c=3%
0
0
0.2
0.4
0.6
Chordwise Station, x/c
0.8
1
Figure 8:
Amplification ratios of TollmienSchlichting waves on biconvex airfoils: present
method (solid lines) compared to COSAL.
9
8
7
6
Present Fit
2500 Hz
2750 Hz
3000 Hz
3250 Hz
3500 Hz
4000 Hz
4500 Hz
5000 Hz
5500 Hz
6000 Hz
7000 Hz
8000 Hz
9000 Hz
10000 Hz
12000 Hz
15000 Hz
17500 Hz
20000 Hz
25000 Hz
30000 Hz
N
5
4
3
2
1
0
0.5
1
1.5
2
2.5
x (m)
Figure 9:
Amplification ratios of TollmienSchlichting waves a biconvex airfoil: present method
compared to LASTRAC.
The fundamental idea is that a table lookup or regression curve fit can be made by parameterizing the
instability growth rate. For the incompressible case,
similarity solutions are used to relate growth rate
to the boundary-layer shape factor and momentum
thickness Reynolds number. Figure 7 is an example of how the rate of increase of n, the logarithmic
amplification ratio, can be approximated by a linear envelope over a range of unstable frequencies.
Since the Tollmien-Schlichting waves don’t amplify
below a certain Reynolds number, that quantity is
also modeled as a function of shape factor. The end
result is a calculation of n that corresponds to the
most unstable frequency at that chordwise station,
however the value of the frequency is not found.
The present approach builds upon Drela’s para-
Figure 10:
Amplification ratios of TollmienSchlichting waves on a flight test airfoil: present
method compared to LASTRAC.
metric formulas for two-dimensional amplification
rates as functions of boundary-layer shape factor, H,
and momentum thickness Reynolds number, Reθ .
He actually used the so-called kinematic versions,
Hk and Reθk , claiming that these formulations are
less sensitive to compressibility effects up to transonic free-stream velocities. The kinematic formulas
are just the incompressible versions of the H and
Reθ equations where the density is not inside the integrand. The argument is that Tollmien-Schlichting
growth rates depend much more on velocity profiles
shapes than the density-weighted compressible parameters.
For the current work, Drela’s formulation is modified in two ways. First, it was originally designed
around his integral boundary layer method, so many
of the empirical formulas for determining Hk and
Reθk can be dispensed with since those quantities are
calculated directly. Second, and more importantly,
his method does not account for the stabilizing effect of a supersonic free stream, so a new parameter
is introduced. The ratio of wall temperature to the
external, inviscid temperature is added as the third
parameter. The current method uses terms with this
new parameter to modify Drela’s original formulas.
It is based on linear stability results on biconvex airfoils at Mach 1.5 to 2.4.
The parametric fits of amplification rate and critical Reynolds number are expressed in the following
manner:
dn
Tw
(5)
= f Reθk , Hk ,
dReθ
Te
Tw
Reθ0 = f Hk ,
.
(6)
Te
As the boundary-layer equations are solved by the
usual downstream marching, Reθk , Hk and Tw /Te
10
are evaluated at each station. The TS amplification
factor can then be determined by integrating Equation 5 from the point at which the flow exceeds the
critical Reynolds number defined in Equation 6 to
the current chordwise station:
Reθ
Z
nts =
on transition prediction. The next higher step in fidelity from the crossflow Reynolds number method
is, as for the Tollmien-Schlichting case, parametric
fits of linear stability results.
crossflow Reynolds number =
dn
dReθ .
dReθ
crossflow shape factor =
Reθ0
Figures 8 to 10 show n-factors calculated with the
present method compared with linear stability results.
Crossflow Instability Calculation
Crossflow transition is more difficult to model than
the streamwise TS modes. The physical phenomena
is predominantly nonlinear and requires, for greatest accuracy, nonlinear PSE solutions and an additional secondary instability code. Presently, only a
few people in the world can accomplish this for subsonic flows, much less on a supersonic aircraft. Consequently, the methods used during conceptual and
preliminary design of aircraft have been very crude
to date.
Prior to the present work, the most widely used
design-oriented crossflow transition criterion was
based solely on the crossflow Reynolds number, defined as:
ρe |wmax |δcf
,
Rcf =
µe
where δcf , the crossflow thickness, is the height
above the surface where the crossflow velocity component is 1/10th of its maximum, wmax . The
method referred to as the “crossflow Reynolds number criterion” in this work is due to Malik et al.24
who introduced a correction for compressibility and
stated that transition is predicted when
γ−1 2
Rcf = 200 1 +
Me .
2
This formula is based on correlations to experiments
and computations on cones and swept wings in supersonic flow.
It unfortuately turns out that the crossflow
Reynolds number approach is often too conservative
for transition prediction on a typical natural laminar
flow wing. The primary problem is an over sensitivity to wing root and tip disturbances and to wing
sweep. Since low wing sweep has an adverse effect on
inviscid drag and sonic boom, it is important to be
as realistic as possible in determining sweep effects
crossflow velocity ratio =
ρew max δcf
µe
ymax
δcf
w max
Ue
δcf
w max
y
max
Figure 11: Crossflow profile and parameters used for
simplified crossflow instability growth model.
The parameters used for the stationary crossflow
fit are summarized in figure 11. These parameters
were suggested by Ray Dagenhart25 and used by him
to successfully model growth rates calculated with
an incompressible linear stability code. Crossflow
n-factors are calculated by modeling the amplification rates and integrating them while simultaneously
solving for the boundary layer flow:
wmax
, Hcf
(7)
αδcf = f Rcf ,
Ue
Rcf0 = f (Me )
Zx
ncf = αdx.
(8)
(9)
x0
Hcf is the crossflow shape factor and α denotes the
spatial amplification rate of the instabilities. The
point at which the stationary crossflow modes become unstable, x0 , is determined with a crossflow
Reynolds number criterion using Malik’s compressibility correction,
γ−1 2
Rcf0 = 43 1 +
Me .
2
The value of 43 works quite well from Mach 1.4 to
2.2 on a number of different airfoils at several sweep
angles.
The previous parametric models for incompressible Tollmien-Schlichting instabilities were developed using boundary-layer similarity solutions. Unfortunately, there are no similarity solutions possible for compressible 3D boundary layers. Instead,
to help identify the influence of each parameter on
11
Figure 12: Crossflow velocity profiles used in “similarity sequence” boundary layers: baseline, adjusted
to increase crossflow shape factor, and adjusted to
decrease crossflow velocity ratio.
Model of stationary crossflow amplification rate at Mach 1.8
4
3.5
Amplification Rate, −α
i
3
2.5
2
1.5
1
0.5
0
−0.5
0
0.2
0.4
0.6
0.8
1
X/C
Figure 13: Stationary crossflow amplification rate
for “similarity sequences” showing LASTRAC data
(symbols) and the parametric model.
Model of stationary crossflow amplification ratios at Mach 1.8
7
6
5
N for Hcf
N for basic
N for cfvr
model for Hcf
model for basic
model for cfvr
N
4
3
2
1
0
0
0.2
0.4
0.6
0.8
1
X/C
Figure 14: Crossflow n-factors calculated with LASTRAC (symbols) and the parametric fit.
Figure 15: Crossflow n-factors calculated with LASTRAC (symbols) and the parametric fit on three
airfoils from Mach 1.4 to 2.2.
growth rates, “similarity sequences” are used. A
boundary layer is constructed by repeating the same
profile (streamwise and crossflow velocities and temperature) at each chordwise x station and
√ just scaling the dimensions in proportion with x. Therefore, the boundary layer is made of a sequence of
profiles with constant crossflow shape factor and velocity ratio. This way each parameter can be individually adjusted to capture its contribution to
the magnitude of instability amplification. It is expected that such an approach is reasonable when
linear stability codes are used on the similarity sequence boundary layers because LST eigenvalue solutions are local to each x station. If PSE were utilized, however, this type of boundary layer may result in erroneous solutions because the conservation
of mass, momentum and energy is violated in the
mean flow.
An example is presented in Figures 12-14. Figure 12 shows the crossflow velocity profiles used
for three similarity sequences: the baseline and
modifications that independently vary the crossflow
shape factor and crossflow velocity ratio. Figure 13
presents the resulting amplification rate of the stationary crossflow modes. The linear stability solver
in LASTRAC26, 27 is used to calculate amplification
rates, and a parametric fit corresponding to Equation 7 is devised to follow this data. In this case, an
analytic funtion is found that approximately fits the
LASTRAC data, but a table-lookup method could
also be used. Finally, Figure 14 shows the resulting
integrated n-factors.
It turns out that the three parameters picked by
Dagenhart are sufficient to model crossflow distur-
12
modeled with a composite amplification ratio,
q
2
2
N ∗ = (nts /ntscrit ) + (ncf /ncfcrit ) .
The present transition criteria is preliminary, but
it is hoped that for the cases that apply to laminar supersonic aircraft, the current approach will
suffice. Nevertheless, correlations with a full three
dimensional en stability code28 should be performed
to validate and improve the streamwise and crossflow calculations used here, as well as the method
used to combine them into a 3D criteria.
Figure 16: Comparison of stationary crossflow amplification between COSAL and the present parametric method on a tapered wing.
bances in supersonic flows. Figure 15 shows surprising accuracy for this example parametric fit that is
based on data at only one Mach number. This observation confirms the assertion by Dagenhart that
crossflow instabilities are much less influenced by
compressibility compared to the TS modes.
Since LASTRAC cannot calculate flows on tapered wings, the present fit is developed with only
infinite-swept wing flows. It can be seen in Figure 16
that the fit compares favorably to COSAL28 results
on a tapered wing with 30 degrees of aft leading edge
sweep and 15 degrees of forward sweep at the trailing
edge (the small-scale test blade flown on an F-1529 ).
3D Transition Criteria
The interference between TS and crossflow instabilities is still being debated. Dagenhart25 briefly
discusses the issue. He cites Pfenninger’s previous
work and says that as long as both TS and crossflow modes are not simultaneously strongly amplified, their interaction can be ignored. Indeed, when
the disturbances are small in amplitude and accurately described by linear theory there would be no
interference. But Dagenhart does state that, “relatively weak oblique TS waves can distort and stretch
crossflow disturbance vortices to produce rapid, resonancelike amplification and transition.”
The current approach to this problem is to maintain simplicity because it is not well understood.
Two methods have been used. In one, the crossflow and TS modes are taken as completely separate,
with their own critical n-factor (9 for TS and 5 for
CF). In the second method, a slight interference is
Example Design Application
As an example application of the boundary layer and
transition estimation methods, a simple wing geometry was optimized to identify the minimum total
drag design and to illustrate the trade-off between
inviscid and viscous drag components. The example
also served to evaluate the feasibility of optimization
tools that include viscous drag and transition modeling for more complex wing design work in the future. By employing a direct optimization approach,
the trade-offs between pressure gradient effects on
TS and CF transition and on inviscid wave drag and
structural constraints may be incorporated in a rational way.
The software used in this exercise consisted of a
number of computer codes including an inviscid Euler solver, boundary layer computations, transition
and skin friction calculations, a response surface generation routine, and a nonlinear search method, integrated into an automated system and executed on an
inexpensive cluster of personal computers. The wing
inviscid (Euler) analysis was completed using Jameson’s FLO107 with a refined grid for accurate drag
and inputs to the boundary layer analysis. A parallel (MPI) version of this code was used to provide
reasonable turnaround times. Next, boundary layer
properties were computed based on the Euler results and the new sweep/taper analysis code. Properties of the boundary layer were computed along
12 arcs across planform. From this data, cross-flow
and TS transition were estimated based on the previously described parametric fit of LASTRAC data
using a combined transition estimate based on an
interaction-free model. Finally the skin friction drag
is estimated using results from the laminar analysis
and a Cebeci-Smith model for the turbulent boundary layer which is added to the inviscid drag from
the Euler code.
13
Figure 17: Section geometry variations based on
simple parameterization.
In this study a response surface / trust region optimization method was employed to avoid potential
problems with non-smooth viscous results. This approach involves generating quadratic models of the
objective function and constraints, then applying
a constrained nonlinear optimization procedure to
identify the optimal result based on the fit. A revised fit is then computed centered around the new
estimated optimum point and the process repeated.
Hundreds of cases were analyzed to obtain multiple
quadratic fits in the course of the optimization, providing detailed boundary layer data in each case for
12 sections over span.
The ten design variables included angle of attack
and twist along with several section properties at
the root and tip stations. A linear variation of properties between root and tip was assumed. Section
geometry variations included changes in thicknessto-chord ratio, location of maximum thickness, a parameter that affected forward and aft curvature, and
a parabolic camber amplitude. Figure 17 illustrates
the variations in thickness distribution represented
with these degrees of freedom. The objective was
to minimize total drag while constraining the overall wing CL and maintaining a constant section moment of inertia at root and tip to maintain torsional
stiffness.
Figures 18- 19 include selected results for this
example wing with a 40 deg leading edge sweep
at Mach 1.6. In Figure 18, the optimized wing’s
amplification factors for TS and cross-flow are depicted. This wing, which yields minimum total drag,
achieves almost full chord laminar flow over much of
the upper surface, while obtaining only 14% laminar
flow on the lower surface in a compromise to maintain section torsional moment of inertia and maintain low inviscid drag.
The optimizer achieves the large upper surface
laminar extent by reducing the curvature on the forward portion of the sections, carefully trading TS
Figure 18: TS and cross-flow amplification factors
on upper surface of optimized wing.
Figure 19: Typical section Cp distribution for optimized wing.
and CF growth. The lower surface is very stable to
TS, but would need some type of control for crossflow stability on the forward part of section. The
Cp distribution, shown in Figure 19 shows the very
linear upper surface gradient over the full span.
This exercise was repeated to define the wing geometry with minimum inviscid drag. Not surprisingly, this calculation, skipping the boundary layer
analysis, led to a wing with essentially biconvex sections and a slightly higher inviscid L/D than the
wing optimized for minimum total drag. When the
viscous analysis was performed on this optimal inviscid wing some laminar flow was still predicted,
but the total (viscous) L/D was estimated to be 8%
lower.
14
Conclusions
A design-oriented aerodynamic analysis was developed to permit direct optimization of supersonic
wings employing laminar flow. This method, based
on a modified sweep/taper concept and fits to linear
stability theory results, permits rapid computation
of boundary layer properties including estimates of
transition and total drag. Although approximate, it
provides important guidance in the advanced design
of laminar wings.
The boundary layer analysis, coupled with an Euler solver and numerical optimization method, was
used to design a wing with minimum total drag at
Mach 1.6. Extensive natural laminar flow was possible on this wing even with relatively high leading
edge sweep of 40 degrees.
Future work will focus on improving transition correlations, more efficiently handling the very
nonlinear optimization problem, and applying the
method to a variety of design problems.
Acknowledgments
The authors would like to thank Jeffrey Viken of
Innovative Aerodynamic Technologies, Inc. for providing some of the en results; Chau-Lyan Chang
at NASA Langley, author of LASTRAC who provided help with his code and general advice on transition prediction; David Lednicer at AMI who generated Euler solutions of the F-15 test wing, and Peter Bradshaw who supplied sweep/taper programs,
documentation and answered many questions. Discussions with Profs. William Saric and Helen Reed
of Arizona State University, Prof. Eli Reshotko,
and Drs. Peter Hartwich and Art Powell of Boeing also have been especially helpful. These studies
were initially and continually inspired by our interaction with Dr. Richard Tracy and ASSET Research
Corporation, interested in the application of laminar flow to efficient supersonic aircraft. Much of
this work has been supported by Directed Technologies, Inc. and their critical role in this research is
gratefully acknowledged.
References
[1] Manning, V. M., and Kroo, I. M., “Multidisciplinary Optimization of a Natural Laminar
Flow Supersonic Aircraft,” AIAA Paper No. 993102, Proceedings of the 17th Applied Aerodynamics Conference, 28 Jun. 1999.
[2] Sturdza, P., Manning, V. M., Kroo, I., Tracy,
R., “Boundary Layer Calculations for Preliminary Design of Wings in Supersonic Flow,”
AIAA Paper No. 99-3104, Proceedings of the
17th Applied Aerodynamics Conference, 28
Jun. 1999.
[3] Sturdza, P., “Design Optimization of Supersonic Natural Laminar Flow Aircraft Including Transition Prediction,” Ph.D. Dissertation,
Stanford University, 2003.
[4] Bradshaw, P., Mizner, G. A., and Unsworth,
K., “Calculation of Compressible Turbulent
Boundary Layers on Straight-Tapered Swept
Wings,” AIAA Journal, Vol. 14, Mar. 1976, pp.
399-400.
[5] Bradshaw, P., Unsworth, K., “An Improved
Fortran Program for the Bradshaw-FerrissAtwell Method of Calculating Turbulent Shear
Layers,” Imperial College Aero Report 74-02,
Feb. 1974 (Revised Nov. 1980).
[6] Kaups, K, and Cebeci, T., “Compressible Laminar Boundary Layers with Suction on Swept
and Tapered Wings,” Journal of Aircraft, Vol.
14, Jul. 1977, pp. 661-667.
[7] Ashill, P. R., and Smith, P. D., “An Integral Method for Calculating the Effects on Turbulent Boundary-Layer Development of Sweep
and Taper,” Aeronautical Journal, Vol. 89, Feb.
1985, pp. 43-54.
[8] Cebeci, T., and Bradshaw, P., Physical and
Computational Aspects of Convective Heat
Transfer, Springer-Verlag, New York, 1984.
[9] Rumsey, C., Sanetrik, M., Biedron, K., Melson, N., and Parlette, E., “Efficiency and Accuracy of Time-Accurate Turbulent Navier-Stokes
Computations,” Computers and Fluids, Vol. 25,
No. 2, 1996, pp. 217-236.
[10] Jameson, A., and Alonso, J., “Automatic Aerodynamic Optimization on Distributed Memory Architectures,” AIAA paper 96-0409, 34th
Aerospace Sciences Meeting and Exhibit, Reno,
NV, Jan. 1996.
[11] Agrawal, Shreekant, and Powell, Arthur G.,
“Supersonic Boundary-Layer Stability Analysis
of an Aircraft Wing,” Journal of Aircraft, Vol.
28, Nov. 1991, pp. 721-727.
15
[12] Goodsell, A. M., “Computational Analysis of
Blunt, Thin Airfoil Sections at Supersonic and
Subsonic Speeds,” Ph.D. Dissertation, Stanford
University, May 2002.
[24] Malik, M. R., Balakumar, P., and Chang, C.L., “Linear Stability of Hypersonic Boundary
Layers,” Paper No. 189, 10th National AeroSpace Plane Symposium, 23-26 Apr. 1991.
[13] White, Frank M., Viscous Fluid Flow, 2nd ed.,
McGraw-Hill, New York, 1991.
[25] Dagenhart, J. R., “Amplified Crossflow Disturbances in the Laminar Boundary Layer on
Swept Wings with Suction,” NASA TP 1902,
Nov. 1981.
[14] Chang, C.-L., Malik, M. R., “Non-Parallel Stability of Compressible Boundary Layers,” AIAA
paper 93-2912, 24th Fluid Dynamics Conference, Orlando, FL, July 1993.
[15] Malik, M. R., Li, F. and Chang, C.-L.,
“Crossflow Disturbances in Three-Dimensional
Boundary Layers: Nonlinear Development,
Wave Interaction and Secondary Instability,” J.
Fluid Mech., Vol. 268, 1994, pp. 1-36.
[16] Haynes, T. S., and Reed, H. L., “Simulation of
Swept-Wing Vortices using Nonlinear Parabolized Stability Equations,” J. Fluid Mech., Vol.
405, 2000, pp. 325-349.
[17] Saric, W., personal communication, Arizona
State University, Tempe AZ, 2001.
[18] Mack, L. M., “Boundary-Layer Linear Stability
Theory,” AGARD Rept. 709, 1984, pp. 3-1 to
3-81.
[19] van Ingen, J. L., “A Suggested Semi-Empirical
Method for the Calculation of the Boundary
Later Transition Region,” Tech. Rep. VTH-74,
Inst. of Tech., Delft, 1956.
[20] Smith, A. M. O. and Gamberoni, N., “Transition, Pressure Gradient, and Stability Theory,”
Douglas Aircraft Report ES-26388, 1956, also
Proc. Ninth Internat. Cong. Appl. Mech., Vol.
4, 1957, pp. 234-244.
[21] Mack, L. M., “Transition Prediction and Linear
Stability Theory,” AGARD Conf. Proceedings
224, Paris, 1977, pp. 1-1 to 1-22.
[22] Drela, M., “Two-Dimensional Transonic Aerodynamic Design and Analysis Using the Euler Equations,” Ph.D. Dissertation, MIT GTL
Rept. No. 187, Feb. 1986.
[23] Gleyzes, G., Cousteix, J., and Bonnet, J. L.,
“A Calculation Method of Leading-Edge Separation Bubbles,” in Numerical and Physical Aspects of Aerodynamic Turbulent Flows II, edited
by T. Cebeci, Springer-Verlag, New York 1984.
[26] Chang, C.-L., “The Langley Stability and Transition Analysis Codes (LASTRAC) – Part I:
LST, Linear & Nonlinear PSE for 2D, Axisymmetric and Infinite Swept Wing Boundary Layers,” NASA TM, in preparation, 2003.
[27] Chang, C.-L., “The Langley Stability and Transition Analysis Codes (LASTRAC): LST, Linear & Nonlinear PSE for 2D, Axisymmetric and
Infinite Swept Wing Boundary Layers,” AIAA
paper 2003-0974, 41st Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 2003.
[28] Malik, M. R., “COSAL—A Black Box Compressible Stability Analysis Code for Transition Prediction in Three-Dimensional Boundary
Layers,” NASA CR-165925, May 1982.
[29] Kroo, I. M., et al., “Natural Laminar Flow
for Quiet and Efficient Supersonic Aircraft,”
AIAA paper 2002-0146, 40th Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan.
2002.