Sep 17, 1996 - Association, The American City Building, Suite 212, Columbia, MD 21044, (410) 730-2656. Work reported ... Supersonic Aircraft Configurations.
NASA-C_-202390
Research
Aerodynamic Supersonic Adjoint James
Aircraft
Reuther,
Technical
Optimization
on Parallel
Juan Jose AIonso,
Mark
Report
at the 6th AIAA/NASA/ISSMO Optimization, AIAA
for Advanced Computer Science NASA Ames Research Center
Configurations
Formulation
RIACS
Presented
Shape
Institute
J. Rimlinger
96.17
Symposium September paper
96-4045
September
via an
Computers
and Antony
Jameson
1996
on Multidisciplinary 1996,
of
Analysis
and
Aerodynamic Supersonic Adjoint James
Shape
Aircraft
Formulation Reuther,
Juan
Jose
Optimization
Configurations on Parallel
AIonso,
Mark
J. Rimlinger
and
of via an
Computers Antony
Jameson
The Research Institute of Advanced Computer Science is operated by Universities Space Research Association, The American City Building, Suite 212, Columbia, MD 21044, (410) 730-2656
Work reported herein was sponsored Space Research Association (USRA).
by NASA
under contract
NAS 2-13721
between
NASA
and the Universities
Aerodynamic Shape Optimization of Supersonic Aircraft Configurations via an Adjoint Formulation on Distributed Memory Parallel Computers J. Reuther* Research Institute for Advanced Computer Science NASA Ames Research Center, MS 227-6 Moffett Field, California
94035,
J. J. Alonso Department
U.S.A.
t
of Mechanical and Aerospace Princeton University
Princeton,
New Jersey
08544,
M. J. Rirnlinger Simco NASA
Ames
Research
U.S.A.
t
Center,
Moffett Field, California
Engineering
MS 227-6
94035,
U.S.A.
A. Jameson; Department
of Mechanical
and Aerospace
Princeton University Princeton, New Jersey 08544,
ABSTRACT
work describes
dynamic shape aircraft design.
the application
optimization The design
the use of both control tributed
memory
the adjoint
U.S.A.
new multiblock ference,
This
Engineering
method process
theory
Control
equations
theory-based
aero-
to the problem of supersonic is greatly accelerated through
and a parallel
computers.
differential
of a control
implementation
theory
is employed
whose solution
on disto derive
allows for the evalu-
ation of design gradient information at a fraction of the computational cost required by previous design methods [13, 12, 44, 381. The resuiting problem
is then implemented
on parallel
distributed
memory
implementation.
we also presented
Furthermore,
preliminary
during
the same con-
results demonstrating
that this
basic methodology could be ported to distributed memory parallel computing architectures [24]. In this paper, our concern will be to demonstrate
that the combined
be used routinely
power
in an industrial
to the case study of the design
of these new technologies
design
environment
of typical
figurations. A particular difficulty propulsion/airframe integration.
can
by applying
supersonic
transport
of this test case is posed
it
conby the
INTRODUCTION
architectures using a domain decomposition approach, an optimized communication schedule, and the MPI (Message Passing Interface)
To realize
Standard
very
is a need not only for accurate
fluid
but also for design methods capable of creating new optimum configurations. Yet, while flow analysis has matured to the extent that
for portability
and efficiency.
rapid aerodynamic
design
dynamics
(CFD).
methods
The final result achieves
based on higher
order computational
In our earlier studies, the serial implementation of this design method [19, 20, 21, 23, 39, 25, 40, 41, 42, 43, 9] was shown to be effective
for the optimization
complex
aircraft
and the Euler equations[39, method
of airfoils,
configurations
was extended
using
wings,
both
25]. In our most recent
to treat complete
wing-bodies,
the potential
aircraft
and
equation
paper, the Euler
configurations
via a
the potential
Navier-Stokes
calculations
plex configurations, to treat moderately Existing
of CFD
to'produce aerodynamic
are routinely
superior
designs,
prediction
carried
there
algorithms,
out over very com-
CFD-based design techniques are just beginning complex three-dimensional configurations.
CFD analysis
methods
have previously
been used to treat
the design problem by coupling them with numerical optimization methods [13, 12, 44, 38]. The essence of these methods, which
*Student
Member
AIAA
incur heavy
?Student
Member
AIAA
/Student
Member
AIAA
optimization procedure is used to extremize a chosen aerodynamic figure of merit which is evaluated by the given CFD code. The con-
§James S. McDonnell AIAA Fellow °Copyright All rights
(_) 1996
reserved
Distinguished
University
by the American
Institute
Professor
of Aerospace
of Aeronautics
Engineering,
and Astronautics,
Inc.
computational
expense,
figuration
is systematically
variables.
Most of these optimization
ation of the gradient
modified
of the cost function
is very
through
simple:
a numerical
user specified
procedures
require
with respect
design
the evalu-
to the specified
design variables. Thesimplest of the
methods
to obtain these neces-
sary gradients is the finite difference method. In this technique, the gradient components are obtained by independently perturbing each design variable with a finite step, calculating of the objective
function
of the differences. tion algorithm
the corresponding
using CFD analysis,
The gradient
to calculate
is used by the numerical
a search
direction
conjugate
gradient, or quasi-Newton
minimum
or maximum
optimiza-
using steepest
techniques.
of the objective
value
and forming the ratio
function
ally expensive
because
perturbations
optimization
the flow must
in every design
when compared verse methods,
variable.
with other traditional
finding the
along
the search
since it permits
sign of supersonic
aircraft
[13].
They applied
OF THE
The
properties
aerodynamic
calculated
Nevertheless,
it is attractive
Cliff,
two
ingredi-
of airframe/nacelle
ADJOINT which
EQUATIONS
define
the cost
function
w, and the physical
which may be represented
I are
location
by the function
of
.T'. That is,
for I = l(w,Jf)
such as infigure
to two-dimensional
and a change
in F results
Hicks and Van Dam, this method
was
in a change
bl = Olr 6w + OIT 0--_ -d-ff
profile
design governed by the potential flow equation. The method was quickly extended to wing design by Hicks and Henne [12]. Later, in the work of Reuther,
in real time.
these
the effects
of the flow field variables,
the boundary,
be repeatedly
any choice of the aerodynamic
the method
almost
combines
including
FORMULATION
of merit. The use of numerical optimization for transonic aerodynamic shape design was pioneered by Hicks, Murman and Vanderplaats
can be optimized in this paper
ents (adjoint formulations and parallel implementations) to produce a robust, accurate, and efficient method that can be used for the de-
functions
is computation-
design strategies
work presented
descent,
After
strategy
The
integration.
direction, the entire process is repeated until the gradient approaches zero and further improvement is not possible. The finite difference-based
point where configurations
y
(1)
in the cost function. The governing equation R and its first variation express the interdependenceof w and .T"within the flow field domain D:
successfully used for the design of supersonic wing-body transport configurations [38]. However all of these cases, which were confined to finite difference
gradients
on serial
computer
architectures,
(2)
were
limited in their geometric complexity simply due to computational expense. For example, the designs presented in [38] were limited to wing-body
configurations.
Yet it is well known that optimum
considerations
into the design problem
outlined
[38l since the required number of mesh points, which more doubles with the inclusion of nacelles, could not be afforded. In the last few years,
alternative
methods
a Lagrange
multiplier
_/J,we have
perfor_I
mance (especially for supersonic configurations) will require highly tuned nacelle/airframe integrations. It was not possible to include nacelle/airframe
Introducing
for obtaining
-
O1T_w
017"
¢,T
__
in
+
15, 33,
partial
area of research.
exhaustive advantages
Choosing
_p to satisfy the adjoint equation
whereby
14] present
a An
of these approaches
the sensitivity
of some objective
is the adjoint
the first term given by
is eliminated,
function
system
(using
the flow equations) very cheaply,
the same
enables
meaning
mature
with respect to an
each
techniques
as perfected
gradient element
the number
is essentially
eliminated as a constraining factor. Moreover, the adjoint tion (and to a lesser extent the accompanying flow solution) not be highly highly-converged difference
converged
to be useful,
flow solutions
gradients
formulation
timization
of a complete
computers to decrease
are crucial
in computational
of the design problem, configuration
The advent
contrast
to the
to accurate
finite
using distributed the turnaround
memory
cost provided the aerodynamic
still remains
of reliable
is a key enabling
and we find that the desired
= Ol OYr
gradient
_Or [OR] _
(4)
The main cost is in solving the adjoint equation (3). In general, the adjoint problem is about as complex as a flow solution. Therefore, when the number
of design variables
is larger than 2, it becomes
pelling to take advantage of the cost differential between solution and the large number of flow field evaluations determine obtained,
the gradient
by finite
_ can be provided to obtain
differences.
to a variety
an improved
Once
com-
one adjoint required to
equation
of numerical
(4) is
optimization
design.
by op-
MULTIBLOCK
FLOW
SOLUTION
a formidable
and efficient
time of these design
(3)
is only a function of F. Since (4) is independent of 6w, the gradient of I with respect to an arbitrary number of design variables can be determined without the need for additional flow field evaluations.
algorithms
In spite of the large decrease
task.
in significant
soluneed
[45].
the adjoint computational
which
for
to be calculated
of design variables
_r
0---_
formulation
arbitrary number of design variables is obtained with the equivalent of only one additional flow calculation. Here, the solution of the adjoint
=
report on the various approaches to the problem and their and disadvantages is given by the first author in [45].
The most promising
,-.
design
29, 16, 35, 30, 28, 34, 47, 31, 36, 37, 49, in this developing
_.T)
than
sensitivities have been developed which greatly reduce the computational cost of optimization. References [1, 2, 3, 5, 4, 7, 8, 6, 32, list of recent works
OR
parallel
technology
calculations
to the
In order dimensional tions,
to extend work
the
the single-block
replaced.
As with
methods
to the treatment flow
solver
the single-block
presented
in our earlier
of complete used solver,
aircraft
three-
configura-
in [21, 40, 421 must the more
general
be
flow
solver must
meet fundamental
requirements
and robust convergence to be employed environment. High accuracy is required provements
in the design
as the accuracy
realized
by the method
of the flow analysis.
is also critical
since
of accuracy,
2ncl
[_=',,_I
Halo
can only be as good
Efficiency
the optimization
efficiency,
in an automated design since the predicted imof the flow solver
of the design
will generally
require the computation of many flow solutions or other solutions of comparable complexity. Finally, robust convergence is also of significant
importance
mization efficiency. noise
since the main
is in obtaining The solutions
in the figure
benefit
of aerodynamic
the last few percentage must be converged
of merit,
the level of realizable
well enough
code written
improvement.
The desirable
achieves
and residual converge was
single-block
by the fourth author
ria. FLO87
ability
smoothing.
to meet
multiblock
these
c]]_gct
the FLO87
DJ-_ti_
ire_ St_T_=il
easy to obtain
solutions
that
in the present
within
the framework
The use of a multiblock
approach
Figure
1 : Convective
and Dissipative Stencils.
work of a block.
approach.
First, the blocks that comprise external file. Then, the double
also currently
under
The general to construct
alternatives
such
as unstructured
mesh solvers
are
investigation.
strategy
in developing
and update
the multiblock
a halo of cells
around
flow solver is
each block
such that
The
strategy
for a complete
each individual grid metrics, For the coarse is repeated
with coarse
ceils adjacent
of adjacent
coarse
with appropriate
boundaries
the halo cell values
flow field data at the appropriate
plish this task, a two-level The
and loading
requirement
times.
halo is constructed
of this double
halo results
To accom-
around from
each
block.
the necessity
of
preserving a complete stencil of calculated fluxes entering and leaving each cell in the entire domain without regard to block boundaries. This ensures
that the conservative
maintained.
Since
are calculated
both the convective cells
are necessary,
level halo for each
of the control
thus requiring
block in the multiblock
The dissipative
fluxes are composed
der differences
corresponding
derivatives
algorithm
fluxes
volumes),
the existence calculation.
of a blend of first and third or-
to terms that mimic
of the flow quantities
is fully
and the dissipative
at the cell faces (boundaries
all six neighboring of a single
flow solution
second and fourth
symmetry,
of convective
and dissipative
each block, some of these cells will lie directly boundary, in which case, the values a different block will be necessary
ceils that are fluxes.
For
next to an interblock
plied. methods
depending
on the kind of boundary
Once the halo configuration for spatial discretization
ceils
in the multigrid
procedure,
grid halo cells defined
grid blocks.
or far field boundaries,
For block
for
blocks.
the process
by the internal
faces
standard
values
of adjacent
cells
that lie on solid,
single-block
techniques
are used to define the halo cells. As an example, consider the simple 4-block grid depicted in Figure 3. The halo cells for block I will be obtained from the internal cells of blocks II, 111, and IV, and from solid
or far field boundary
other
blocks.
aggregating cell process. all coarse
Coarse
techniques
for the faces not adjacent
grids are computed
in the usual
groups of eight cells and then repeating Once the halo configuration is complete
grids, the flow solution
fashion,
to by
the above halo for the fine and
commences.
The system of equations solved here as well as the solution strategy follows that presented in many earlier works [26, 18, 17]. The threedimensional
Euler equations
may be written
as
Ow Of, 0--7 +_=0 where
it is convenient
locity components defined as
to denote
inD,
(5)
the Cartesian
coordinates
and ve-
by Jcl, x2, z3 and ul, u2, u3, and w and f,
are
of the flow variables residing in to calculate the convective and
dissipative fluxes. Halo cells on the external boundary of the entire computational domain are constructed and updated by extrapolation and reflection,
as follows:
[261. For the third order differences
each cell within a block, Figure I shows the neighboring for the calculation
proceeds
into halo cell locations
from the interior
grids required
at the boundary faces of each cell for all blocks, the presence of the twelve neighboring cells (two adjacent to each face) is required. For required
flow solution
the flow field mesh are read from an halo configuration is established, for
block, by inserting etc., taken
the flow solution inside each block is transparent to the block boundaries. This task requires establishing the size and location of halo to block
Discretization
is a first
step towards the treatment of more complex configurations. However, the multiblock strategy presented here is not the only viable Other
[{_Tci i
met all of the above crite-
The challenge
swict requirements
flow solver.
also
with the aid of multigridding
It is normally accuracy.
to com-
as a check
applications,
readily
fast convergence
to machine
that the
say drag at a fixed lift, is well below
pare adjoint-based gradients with finite differencing requires highly-converged flow solutions. In our three-dimensional
opti-
points in improved
condition
is set up for each block, and time integration
ap-
standard
(including
ar-
pul W
=
pu2
,
f,
=
pU3
pu,ul
+ p6,1
pUtU2
"_-p6,2
pit,U3
-_-
oE where 6,jis the Kronecker
pu, H delta.
tificial dissipation, implicit residual averaging, and multigridding) are employed to compute the flow solution within each individual p=(7-
(6)
p613
1)p
Also,
{ E-_ 1
(7)
and
pH = pE + p, where 7 is the ratio of the specific heats. to coordinates
Consider a transformation
j'
to eliminate the need to solve large lxidiagonal systems spanning blocks, which would incur a penalty in communication costs may not even be defined.
_1, _2, _s where
/¢"=/ae, Introduce
(8)
converged significant
J=det(K),
scaled contravariant
velocity
Ui
_-
Q'3
L°_' J.
K,;' components
solution, reduction
This
change
has no effect
and in the present application in the rate of convergence.
the and
on the final
has not led to any
THE ADJOINT FORMULATION FOR THE EULER EQUATIONS
as The application
Uj
of control
theory
to aerodynamic
design
problems
is illustrated by treating the case of three-dimensional design, using the Euler equations discussed above as the mathematical model for
where Q= The Euler equations
Jtf
-1 "
can now be written
compressible flow. In our previous work, the illustrative problem most often used specified the cost function as a measure of the differ-
as
ence between OW
OF,
0---t-+ _
= 0
in D,
(9)
with
configurations,
pul
pU, ul + Q,lp
pu2
, F, = Q'sfJ
=
pU, u2 + Qi2p
•
pl13
,
J
pE
U n=0
where
boundary
for incoming solution.
a solid
faces,
waves,
while
The time integration strategy
[26].
The
balance,
updating
BF,
freestream
outgoing
scheme
solution
follows proceeds
the flow variables,
that
is normal
to exist.
conditions
waves
distributions
problem
for supersonic
design:
---_
C)' D
=
CA cos c_ + C,v sin
=
1
[[
SrefJJ.s
CP (Sxcosa
+ Sysinot)
At the
are specified
are determined
by the
plane of the face in question.
Note
that the integral
in the final
expression above is carried out over all solid boundary faces. The design problem is now treated as a control problem where the control function
is the geometry
I, subject
shape,
to the constraints
which
is to be chosen
to minimize
defined by the flow equations
variation in the shape will cause a variation consequently a variation in the cost function
by performing and smoothing
bl =
the cell flux the residuals,
_CAcosc_
+6CNSino_
{--CA
+
[ OCA cos a + -_OCN sin a } 6a t-'_--a
each stage of the integration
configuration
to be used, without
permits
process.
standard
modification,
time stepping,
The only implementation single-block
difference
and multigrid
where
single-block
for the computation
between
of the implicit solution
strategy,
convergence
the integration residual
_
of
6CA
and _CN
are variations
due to changes
in the design
parameters with e_fixed. To treat the interesting problem of practical design, drag must be minimized at a fixed lift coefficient. Thus an additional
constraint
is given by
acceleration.
strategies
averaging
a tridiagonal
sin o_ + Ct¢ cos_}
The
the flow field within each individual block. This includes the singleblock subroutines for convective and dissipative flux discretization, multistage
A
that used in the single-block +
subroutines
(5-10).
6p in the pressure and
at
use of the double-halo
at a fixed
d_ld_2,
each stage of the time stepping scheme and each level of the multigrid cycle. The main difference in the integration strategy is the need to loop over all blocks during
that will
on
(I1)
is assumed
of pressure
where Sx and S v are projected surface areas, Sre f is the reference area, and d(1 and d(2 are the two coordinate indices that define the
on the direction
surface
the specification
drag minimization
I
conditions
on allBs
0 is 1, 2, or 3 depending
far field
In
salient lift.
notation applies to each block as the steady-state solution to
equation (9) in all blocks subject to the flow tangency all solid boundary faces of all blocks:
Bs
distribution.
is a considerably more chalhere will focus on the more
pU, u3 + QI3p pU, H
For the multiblock case, the above in turn. The flow is thus determined
to the face
pressure
determine near optimum performance lenging problem. Thus the development
(1o)
where
and some desired
optimum performance is well understood by most aerodynamicists. However, for the case of supersonic design of three-dimensional
p
W = J _
the current
the case of transonic flows over conventional commercial transport wings this aerodynamic figure of merit proves to be very effective since the tailoring of these pressure distributions to achieve close to
technique.
system
8CL = 0,
is in the In the
of equations
which gives
is
$6'N COSa -- 6Ca sin a
set up and solved using flow information from the entire grid. Thus, each residual is replaced by an average of itself and the residuals of
+
the entire grid. In the multiblock strategy, the support residual smoothing is reduced to the extent of each
+
for the implicit block, in order
{ --C'N sin c_ -- CA cosoQ OCN cos. -Tff-.
6c_
OC,_ O. sin c_} 6c_ = 0
Combining
these two expressions
to eliminate
_,
gives
Suppose
now that
O is the steady-state
solution
of the adjoint
equation 61
= 6Ca cos,
8¢
+ _C_v sin,
&/, -C'[_=0
+_ { _c_ cos. - _c. sin. }, (12)
At internal jacent
block boundaries,
blocks.
conditions
where l] is given by
depends
on whether
tional
domain
_1-5
on w through the equation of state (7-8),
tion 6p can be determined
is used. the variations
in the mapping
the varia-
from the variation 8w. If a fixed computa-
derivatives.
in the shape result in variations
Define the Jacobian
Of, A, = -ff-ww'
ma_ces
If, however,
C, = QoAj.
(13)
for 6w in the steady state becomes 0
on all
the flow is subsonic
outgoing
waves
(17). For supersonic domain
Now, multiplying by a vector co-state variable q,, assuming the result is differentiable, and integrating by parts over the entire domain,
analysis.
in the surface
at any particular
For a given
over which
local conditions. adjoint
t
_l
J
fi, are components
The variation
JB (r"er6F')
of a unit vector
in the cost function
normal
&5,
(14)
to the boundary.
can also be expressed
6p after (12) and (14) are summed
is fairly close,
0 for incoming
physical
geometry,
point.
surface
changes
waves,
It is noted
conditions
at the
intuition
as well as
say a wing,
a change
P, will
incur changes
in the
in terms of
the flow condition
of this reverse
The contribution
by changes
affect
at a
also form a cone that would point of the Mach cone, depending on
It is the solution
represents.
influenced __
=
by the solution.
from
given point. This region would exactly in the opposite direction
--lit.
BE.
flows, the choice of boundary
can be developed
mathematical
the region
where
' a,,-" ._ ' &%
is very far from the
and the boundary
are determined
boundary
or supersonic.
flow field and hence the performance in the region defined by the Mach cone originating at "P. Similarly. it is possible to determine
0,,-7(_F,) = 0,
where
=0
from the ad-
that the waves in the adjoint problem propagate in the opposite direction to those in the flow problem due to the transpose in equation outer
Then the equation
the flow is subsonic
then far field faces may be set by _,1-5 while
cancel
of the adjoint
flow, so long as the outer domain of interest, we may set
(Or + _a_ cos. - °cAo,, sin .) ' Since p depends
(17)
the face integrals
At the far field the choice
For subsonic configuration
= (-C_+_cos,+_sin,)
inD.
to the surface
problem
to. say, drag at a given at all points
within
that the point is
the reverse
cone. The correct supersonic far field boundary conditions for the adjoint equation that are consistent with this reversed character are: _],l-s
=
0
_&-s
at the exit; extrapolated
from the interior
at the inflow boundary.
to give, Then if the coordinate transformation is such that 6 (JK -1) is negligible in the far field, the last integral in (15) reduces to
_'_M_Sref
s +f2 (Sy cos c_ - S_ sin a.) }
-- /fBs Thus
+
I Sref
ffB
s
C T 6 F'7 d_ld_:.
(18)
d_ld_2
by letting
_/, satisfy
the boundary
conditions,
Cv{(6S'_c°sa+6Susin_) (02Q,71
+ q_3Q,72 +
_/,4Qrt3)
:
Q on all Bs,
(19)
+f2 (6Sy cos _ - 6Sx sin a) } d(, d(2 where
- JDI ( On the solid surfaces (1 1) that
Bs,
,< +
(,,) Q
hi =
fi2 = 0. It follows
1 , _TM_Sref
-
from equation
,
,
+ S_sin.)
by parts again that
5 (Q,7,)
0 {0 Qo26p
cos,
cos.- s=sin.)}, we find after integrating
Q,7, _p
{(S.
. + p
6 (Qn:)
Q,738p
_5(Q,73)
0
0
S
on any
Bs.
+ (16)
which is independent
(20) of w.
MULTIBLOCK ABLES
MESH
VARIATIONS
AND
DESIGN
VARI-
Thus in order to use analytic
61 in equation (20), the variation
terms must be obtained to use finite differences avoids
gradient,
but it unfortunately
the use of multiple
used repeatedly.
in the metric
in each block. One way to accomplish to calculate the necessary information.
approach
flow solutions
still requires
The number
to determine
the mesh
of mesh generations
this is This
generator
required
the to be
is propor-
tional to the number of design variables. The inherent difficulty the approach is two-fold. First, for complicated three-dimensional configurations,
elliptic
must often be solved These
or hyperbolic
iteratively
iterative
partial
differential
equations
in order to obtain acceptably
mesh generation
procedures
in
ods for obtaining
metric
variations
smooth
with an iterative
in a way that resembles
make use of the relative
a trivial task. In fact no method
currently
In our earlier
works
automatic
140, 39, 25, 19, 20, 21], two methods
mesh points
the block.
interior
been explored which avoid these difficulties. In the first method, a completely analytic mapping procedure was used for the mesh This technique
is not only fully automatic
and results
in
smooth consistent meshes, but it also allows for complete elimination of finite difference information for the mesh metric terms. Since function
fully determines
the entire
surface shape, this analytic relationship finite
to obtain
step.
the required
An analytic
information
mapping
mesh based on the
may be directly method
differentiated
point
to be highly
effective
without
considering
a
requires
the geometry
[19, 20, 211.
The second method that we have explored is the use of an analytic mesh perturbation technique. In this approach, a high quality mesh appropriate procedure
for the flow solver
is first generated
prior to the start of the design.
by any available
In examples
is shifted
corner
point.
The second
to produce
stage corrects
a very simple are created
algebraic
by moving
mesh
This permitted
perturbation
algorithm.
say the the use of
New meshes
all the mesh points on an index line projecting
from the surface by an amount which is attenuated as the arc length from the surface increases. If the outer boundary of the grid domain is held constant
the modification
to the grid has the form
= ,:,d + so,d
These
perturbations
the volume
the perturbations
resulting
are then also incorporated
Finally
with both corner
for, the third stage checks
into the
and edge point
the perturbation
motion
of each point
in all six faces relative to the position of the corresponding point in the stage 2 block. If the perturbation of the domain involves a simple translation of all boundary points, the relative from stages 2 and 3 will be zero and all the perturbation accomplished
by stage
1.
If, however,
face points
changes will be
are perturbed
relative to the reference block, then these changes are propagated to the interior points through relative arc length-based perturbations along
projecting
index
lines.
In general
all 3 stages
is to use an initial
mesh
are required.
with good
quality
attributes as a starting point, and then systematically perturb this mesh in a manner such that the original grid quality is maintained, without
the need for expensive
Since our current
elliptic
smoothing.
flow solver and design algorithm
to-point match between perturbed by WARP3D,
assume
a point-
blocks, each block may be independently provided that perturbed surfaces are treated
continuously across block boundaries. a new mesh is given by the following
The entire method algorithm.
that touch
in an adjacent 3. All inactive
faces
4. WARP3D
by the design
of creating
variables
(active
face, either in the same block or
are implicitly
that either
perturbed
include
edge or abut to an active face quasi-3D form of WARP3D.
or implicitly points.
- x:',
an active
block,
by (21).
an implicitly
are implicitly
perturbed
perturbed
by a
is used on each block that has one or more explicitly perturbed
faces to determine
the adjusted
interior
(21) Note that much of the mesh, especially
where z i represents
to one
to be shown
that only one surface,
a design case.
proportional
volume mesh points through a weighting scheme that is proportional to the relationship of an individual edge point motion and the volume
2. All edges
during
points
of each corner point,
from the first stage by determining the distance each of the 12 edges of the stage 1 block needs to be moved to attain the desired edge
all subsequent meshes which are developed by analytical perturbations. In the method that was previously developed for wing-body wing, was perturbed
block that
along the index lines away from that
1. All faces that are directly affected faces) are explicitly perturbed.
it had been assumed
an interim
of the 8 corner
to the motion
later, these meshes were created using the Gridgen software developed by Pointwise Inc.[46[. This initial mesh becomes the basis for
configurations
(TFI)
in the initial mesh.
by a displacement
distance
The idea of WARP3D
topology to be built directly into the formulation, and only works for simple configurations. Nevertheless, within these limitations it has proven
Corresponding
minus the normalized
accounted
have
interpolation
point distributions
by the new locations
entirely
point in question.
is by no means
interior
is determined
generation
exists that allows this to be accomplished as a completely process for complex three-dimensional configurations.
transfinite
stage shifts the internal
locations.
in order
for the
The WARP3D algorithm has been modified from that presented in reference [431 and is now a three stage procedure. The first
mesh generator leads to computational costs which strongly hinge on the number of design variables, despite the use of an adjoint solver to eliminate the flow variable variations. Second, multiblock mesh
the mapping
perturbations
1481. Unlike TFI, where there is no prior knowledge of the interior mesh, the perturbation algorithm developed here CO/ARP3D) does
each
are often com-
in combination
be modified
defining
putationally expensive. In the worst case they approach the cost of the flow solution process. Thus the use of finite difference meth-
generation.
mesh
treatment of the more general problem where multiple faces of a given block may be simultaneously deformed, equation (21) had to
In order to construct
meshes.
move.
grid points, x.., represents
the surface
not require
mesh perturbations
away
from the surfaces,
and thus may remain
will
fixed through
grid points and S represents the arc length along the radial mesh line measured from the outer domain, normalized so that ,5 = 1 at the
the entire design process. Close to the surfaces, many blocks will either contain an active face or touch a block which contains an
inner
active face, either by an edge or by a corner. As the design variations affect the active faces, the above scheme ensures that the entire mesh
surface.
case where
Unfortunately multiple
faces
this simple sharing
logic
common
breaks
edges
down
in the
are allowed
to
willremain attached along blockboundaries. Added complexity is needed edges
to accomplish and comers
master
step (2) since the connectivity
must be indicated
edges and master
comers
somehow.
Since
this mesh perturbation
Currently,
is determined
step. During the design calculation, these transferred to all connected edges and comers algorithm
of the various a list of
as a preprocessing
since
the current
the mesh perturbation
paper was significantly
algorithm
that was
more complex,
used in
calculate 6Q,j instead of deriving the exact analytical relationships. Even in cases with hundreds of design variables, the computational re-evaluating/_Q,i
for all necessary
blocks
is
still insignificant compared with the cost of a single flow solution. The conclusion is that the analytical mesh perturbation algorithm, WARP3D,
unlike an elliptical
mesh generation
to the extent that the cost of remeshing
method,
is efficient
to these design
variables
can be neglected.
6I
where
61 is calculated
perturbed
by (20)
8hi'
(22)
and each
by a finite step. Therefore,
and on the appropriate
DOMAIN TATION
DECOMPOSITION
the design
code are: a domain
Program
Multiple
Interface)
Library
Data)
term
G, a basis space of
to homogeneous
reliability.
which were initially
were chosen
as a
simply for their ease of implementation The form of the Hicks-Henne
proposed
in Reference
by domain
[(
=
sin
functions
necessary
The
networks
,
0 < z < 1
(23)
Since
of
choice
Passing
of message
obtaining
the desired
entails the treatment
of four
No attempt is made algorithm. It is thus
when
compared
of the step sizes and algorithm is com-
with the other
elements
the design.
the flow and adjoint
equations
numerical
are to be solved
techniques,
used for the flow equations
using ex-
the same parallelizaapply to the solution
of
the adjoint equations. Therefore, all details of the parallel implementation corresponding to these first two parts of the program will with reference
tain mesh
only
the mesh perturbation basis,
consistency
decomposition
been used for the modification
of
of workstations.
sections,
decomposition
insignificant during
on a block-by-block
a-z_-_(_
which have traditionally
and the MPI (Message
assumed in this construct that the determination search directions provided by the optimization
be explained
[12] is given by:
Here tl locates the maximum of the bump in the range 0 _< z _< 1 at z = tl, since the maximum occurs when z '_ = 1 where 3' = log ½/ log tl, or a log t_ = log _-. The parameter t2 controls the width of the bump. To generalize the application of the Hicks-Henne functions,
the parallelization
separate parts: the solution of the flow equations, the solution of the adjoint equations, the calculation of the mesh perturbations, and the
more, since b(z)
IMPLEMEN-
model, a SPMD (Single
passing.
and heterogeneous
tion techniques
and their proven
strategy,
As one can see from the previous
of the design space.
In this work, design variables
to be
passing library was determined by the requirement that the resulting code be portable to different parallel computing platforms as well as
putationally
b, is independently
to construct
such
are allowed
PARALLEL
decomposition
for message
actly the same efficient
functions
AND
The main strategies that are used to accomplish
independentperturbation functions hi, i = !, 2,..., n (n = number of design variables) must be chosen to allow for the needed freedom set of Hicks-Henne
faces in both blocks
while both surfaces
calculation of the gradient integral formulas. to parallelize the quasi-Newton optimization
as
_(b,)-
and then with
does not change
parallelization
It remains to choose a set of design variables which smoothly modifies the original shape, say b,. The gradient can then be defined with respect
weights
that thickness modified.
of a wing is to be preserved
halves of the wing are in separate blocks, need to be applied at the proper locations
and it was dis-
covered that the computational cost of repeatedly using the block perturbation algorithm was minimal, finite differences were used to
cost of repeatedly
if the thickness
it is possible
to work out the analytical variations in the metric terms required for equation (20). This approach was followed in reference [40]. However
For example,
the upper and lower the design variables the proper
lists are updated and as the mesh is moved.
is analytic
face.
algorithm
WARP3D
the communication
can also be addressed
strategies
The subdomains
to the flow equations.
Furtheralso works
necessary
to main-
by the same domain
that are used for the state and costate
of the flow solution
resident
fields.
on each processor
are
divided along the block boundaries such that one or more complete blocks are allocated to each processor. This has the natural consequence through
that the communication between subdomains the same halo cells that were described earlier
is performed for the multi-
airfoil sections where the z in equation (23) refers to the chordwise position, they are applied directly to a parameterized (fi, _) surface
block computations, only now, the interior cells corresponding to given halo cells might reside in a different processor. An alternative
which
to this would be to partition the complete problem along the three coordinate directions for each of the blocks in the mesh. Since the
may be composed
The parameterization study,
the wing
perturbed
of one or more
faces in different
blocks.
may be accomplished
in many ways.
For this
is designed
onto
a plane
Thus the Hicks-Henne in either
and shape
by projecting normalizing functions
the fi, b, or both directions.
structured
so that these variables
the parametric
surface.
all surface
may be applied The design
may be applied
Alternatives
points
by the planform
are provided
to be
outline.
as functions
code is further to any subset such
of
that these
sizes of the blocks can be quite small, severely
limit the number
the flow and adjoint
responsible
case in which
may be prescribed
to a highly
versatile
to have at least
faces in the multiblock
domain
at the input
may be perturbed constraints,
by the design
each design
level,
one or more variable
leading
variables.
To enforce
may be activated
geometric
on more than one
The underlying
for the computations
direction as opposed All of these options
code in which
solutions.
levels
would
that could
be used in
assumption
is the fact
that there always will be more blocks in the mesh than processors available. If this is the case, every processor in the domain would be
variables may be linearly lofted in the second to Hicks-Henne functions in both directions. design
this further partitioning
of multigrid
blocks
there
can be adequately
the advantage the parallel
as many
split during blocks
that the number implementation
serial version,
inside
one or more blocks.
are more processors
than blocks
a pre-processing
as processors.
of multigrid
the
step in order
This approach
has
levels that can be used in
of the code is always
and is only limited
In the
available,
the same as in the
by the size of the smallest
block in
themesh. Oneadvantage oftheassignment of complete processor
is that the number
an arbitrary
integer.
of processors
The drawback
1. Decompose
blocks to a given in the calculation can be
of this approach
exact load balancing that can be achieved partitioning. All blocks in the calculation
by coordinate direction can have different sizes, will be This, in
turn, will imply that some of the processors
will be waiting until the
processor
has completed
of cells
with the result being that the parallel approach
that we have followed
is to assign number
to each
processor,
performance
of the processor 2. Solve
will suffer.
problem
in a pre-processing
a certain
step,
,
The
of cells in each processor
based on the total number
the size of each
communicated are the subject
block,
play an important of current study.
Now, within each processor, to communicate
and the size of the buffers role in proper
there will be several blocks.
to pack all data that needed
from one processor
to another
the number
regardless
array.
of and
In order
large arrays that need to be sent to all other processors. Similarly, another pointer list for the locations of the data to be received is also setup. At the time of communication exchanges each processor communicates all the information for the blocks that it contains to that
implemented using MPI constructs. of volume
calculated
for the gradient
communication
(see [22[)
is
include
a
that can also be
of the integrations
cation at all.
Once a single
the individual global
is such that.
the need for any communi-
processors
MPI reduce
into the total gradients,
merical optimization design iteration.
without
algorithm
have completed
operation
their
can combine
which can then be passed to develop
equation nonzero
• Determine
with an
locations
via
faces.
6Q,,j,
within
those
by finite differencing.
(20) to obtain
61 for those blocks that
6Q,.:.
the gradient
component
by equation
(22).
the search
direction
and perform
a line search.
has not been reached.
the
to the nu-
the shape change
at each
The basic method
here builds on that used in reference
proper
to treat multiblock
extensions
ment the method,
equation
ALGORITHM
(SYN87-
domains.
(17) and boundary
140] with the
In order
to imple-
condition
(19) must
be discretized on the multiblock domain. In the current implementation, a cell centered, central difference stencil that mimics the flux balancing used for the flow solution is used. Since this choice of discretization
differs from the one obtained
tion Jacobian matrix system, the gradients exactly discrete
ff the discrete
were actually transposed obtained by the present
flow equa-
to form the adjoint method will not be
equal to the gradients calculated by finite differencing the flow solutions. However, as the mesh is refined these differ-
used for the flow solution.
Therefore,
flow solver,
uses an explicit
multistage
accelerated
by residual
communication
smoothing
is again handled
for the full transfer
the adjoint
and multigridding.
through
of information
solver,
Runge-Kutta-like
across
like the algorithm
Intra-block
a double halo which allows boundaries
except
for the
stencil of support for the implicit residual smoothing. Step (3) in the above procedure is the portion of the method is still treated
by finite differences.
incur only a trivial
computational
that
Fortunately,
all of these steps
cost compared
with even a single
flow analysis time step. It is therefore possible, without significant penalty, to leave this in finite difference form even for cases where of design
variables
hanced
by Kennelly
are used.
uses the quasi-Newton algorithm, Murray and Pitfield [10] and en-
[27], to calculate
the search
direction.
It is an
unconstrained optimization algorithm that uses Broyden-FletcherGoldfarb-Shanno (BFGS) updates to the Cholesky factored Hessian
With all the necessary based design,
DESIGN
point
all the delta metric terms,
The present implementation QNMDIE developed by Gill, MULTIBLOCK
mesh
variable.
an edge
blocks with perturbed
that were perturbed
Integrate
many hundreds
COMPLETE MB
by thedesign
ences should vanish. Continuing, the adjoint system so discretized is solved on the multiblock domain in a fashion identical to that
send and receive
of the cost function
integrals
The nature
can be calculated
calculation,
The
decomposition point of view, the problem is nearly parallel, since the partial values of the integrals in
each domain
results
it.
(non-blocking)
and surface
in parallel.
from a domain embarrassingly
need to receive
asynchronous
The final formulas number
•
internal
6. Return to (2) if minimum
Within
adjoint variables,
are stored on a large one-dimensional
processors
blocks
new
for those
• Calculate
5. Calculate
to accomplish this type of communication, during the pre-processing step, each processor compiles a pointer list with all the enmes in these
those
the
all faces affected
to
to be communicated
in one single message,
the data for the flow variables,
grid locations
perturb all faces that share perturbed design variable.
Implicitly explicitly
WARP3D
blocks that need
by a finite step according
The data for these
of blocks that resided in each of the processors.
each processor,
to the boundary condi-
•
• Obtain
and
neighboring blocks can reside in a different processor, and therefore, communication is necessary. In order to minimize communication cost, it was decided
(5-10).
(17) subject
Perturb the design variable (23, etc.).
contain
with their neighboring
equations
equations
perturb
to be
load balancing,
the flow field governing
• Explicitly
of cells in
each processor is only an approximation to the optimal solution of the problem. Other variables such as the number of blocks, the block locations,
numberof
for the communication
4. Por each of the n design variables repeat the following:
as possible to the exact share for perfect load balancing. this approach yields quite good load balancing [24]. One
must note that load balancing
lists of pointers
halo cells.
3. Solve the adjoint tion (19).
its work
to solve the load balancing
of blocks such that the total number
is as close In practice,
and create
is the loss of the
and consequently, it is very likely that different processors assigned a different total number of cells in the calculation. with the largest number
mesh into an appropriate
the multiblock
processors,
components
it is now possible
defined to outline
for the multiblock the complete
adjoint-
procedure:
matrix. A complete treatment of the quasi-Newton and other optimization strategies is given by Gill, Murray and Wright [11].
Mu_b4ock
PARALLEL
EFFICIENCY
For problems
with a low task granularity
3°
nately, base
convergence
domain.
Thus,
techniques mesh
acceleration
their success
current
could
sizes.
complex
on global
However,
geometries,
allel efficiency
meshes
used in viscous
becomes
of the design
parallel
method
redesign
used containing
72 blocks
plementation
is taken
In addition,
solutions.
Finally,
in the performance
calculations.
moderate
as this, excellent
size such
up to 16 processors. shows
a breaking
load balancing
the complete
serial
speed-ups
Just as important
a 72 block
of
calculation good
into 32 processors.
shows a spread in the load with the fact that the serial with short exe-
accounts
completed
in under
curves
Two
for the decrease
35 minutes
of the method
complete
design
is the true
iterations
using 32 processors.
were
The rapid re-
sponse of the design method clearly allows the possibility of investigating a large number of configurations, and increases the probability of arriving
at a truly optimum
design.
Timings
2: Parallel
break point.
for the 1,500,000-cell
HSCT mesh used in the Results sections can be readily extrapolated from the ratio of the numbers of cells in the two meshes.
The symmetric
initial airfoils
NUMERICAL
TESTS
To demonstrate
the utility of the design
algorithm.
AND RESULTS
will
The baseline
be used
over the span up to the leading ing edge break
edge break.
where the wing sweep
method,
transport
a supersonic
trans-
for the optimization configuration
depicted
in Figures 4----6 was sized to accommodate 300 passengers with a gross take-off weight of 750,000 lbs. The supersonic cruise point is Mach
2.2 with a CL of 0.105.
As can be seen in Figure
planform
has a break in the leading
edge sweep.
edge sweep is 68.5 degrees the Mach penalty. range
while the outboard
angle at M = 2.2 is 63 degrees
edge bluntness
may be used inboard
Airfoils
which
from 4% thick
use blunt
30%
with the pur-
and 80% chord
Outboard
of the lead-
is ahead of the Mach
sharp leading
edge was used to avoid undue wave drag.
were chosen
to be symmetric,
biconvex
shapes
cone,
a
The airfoils
modified
to have a
region of constant thickness over the mid-chord. Figure 5 shows that the four-engine configuration features axisymmelaic nacelles tucked close to the wing
lower
drag by minimizing
surface.
This layout
the exposed
because
diverter
of the channel
reduced
However,
flows occurring
region of the diverter,
wing, and nacelle
leading
of the diverters
edge heights
favors
area.
wave
it may be
in the juncture
at the wing trailing edge.
are determined
The
by the boundary
layer local displacement thickness such that entrainment of boundary layer flow into the engines is avoided. Since the distance from the wing leading edge to the diverter leading edge two nacelle, this causes a corresponding diverter In this design problem, ward
we have chosen
task of recambering
the wing
is different for the height difference.
to attempt
the straightfor-
so as to minimize
the total con-
fixed CL and Mach number. than allowing both thickness
and camber to vary while maintaining spar constraints and adequate fuel volume. However, since the complex nacelle geometry and re-
as a testbed
supersonic
Design
were chosen
thick spars at roughly
figuration wave drag while maintaining This represents a much simpler problem
port configuration
i 30
Speed-ups for the AdJoint Method
pose of accommodating
problematic
time.
i 25
was sev-
using
of achieving
mesh
of processors),
as the speed-up
execution
i i 15 20 Number Of Procossonl
execution
are obtained
for the 32 processor
Figure
im-
efficiency.
wall-clock
i 10
was
calculations) which include
play a large role for problems
cution times (large numbers in parallel
topology
I/O time was also included
The optimum solution for this problem balancing of almost 10%. This, together parts of the algorithm
i 5
jet con figuration.
efficient
trend due to the impossibility
by partitioning
•
can
As one can see, for a problem
The speed-up
Q.
The par-
speed-ups below are comThe execution time of the
of the design problem (for parallel performance taken to be that of 2 complete design iterations, eral flow and adjoint
.i/
in Figure 2 for the
as that of the most
of the method.
::
sizes with a total of 750,000
It must be noted that the parallel using true execution wall-times. job
i
25
and the
schedules
a wing-body-naceUe
of varying
C_rw
'1'15
flows and
performance.
is presented
.................................................................. --Ideal i
10
of a typical business
of the aircraft,
in traditional
low enough,
high. of the communication
Sp4NKlup
i+
smoothing
for larger
In the modeling
one processor
residual
performance
with excellent
case of the automated
cells. puted
and implicit
- Parallol
in the 19g0s
in the computational
parallel
the granularity
method
developed
hinder
parallel performance is quite A careful implementation yield a design
techniques communication
multigrid
possibly
of floating point operations can be obtained. Unfortu-
Code
i i i "............. '.....
(ratio of the number o f bytes
received by a processor to the number it performs), large parallel efficiencies
Demgn
6, the leading
is 49.5 degrees.
Since
it is clear that some leading
without leading
at the root to 2.5%
The inboard
a significant edges thick
were
wave drag created
at the leading
that edge
sulting
meshes
on the wing lower
with the wing motion, nificant
challenge
benefits
could
surface
must be moved
even a camber-only
for demonstration
also be attained
design
purposes.
Further
not only by allowing
in unison
represents
a sig-
optimization thickness
vari-
ations in the wing optimization but also by performing both body optimization and nacelle orientation and position optimization. Such studies will be made in future work. Returning shown
to the wing
camber
in Figure 4 was created
putational
cells including
forces a larger number cial from the standpoint
design
at hand,
with 180 blocks
halos. of blocks
While
the mesh
system
and 1,533,440
the point-to-point
to be used, this is actually
of parallel
load balancing.
com-
topology
Ninety
benefidesign
variables
of the Hicks-Henne
and chordwise
direction.
type are lofted in both the spanwise
filling the space
These are spread over the entire wing such
that complete freedom is allowed with the exception of the wing root intersection and the wing trailing edge. The former was held fixed
as the complete
since the capability of reintersecting different geometry elements has yet to be implemented. The latter was held fixed since it was determined
that freedom
in the wailing
waves in the spanwise direction violate possible flap constraints.
edge would
likely
muitiblock
tions that were imposed
produce
along the trailing edge which would Figure 6 shows the solution on the
on the wing.
in the optimization
This allowed
without
the nacelles
approximations.
Due to a shortage
of both resources
to be
computational
in Figure 4
upper surface
pressure
and the Cps
regions
around
This creates wave drag.
since this approach
have been tuned
to optimal
used
with
using both the analytic
ume implementations
on the wing lower
grant
greatly
number
from the generous
AFOSR-91-0391, USRA through Ames Research
support
ARPA
under
of the grant
RIACS, the High Speed Center, and IBM. Con-
Saunders
and M. E. Eleshaky.
of Sterling
Aerodynamic
[21 O. Baysal and M. E. Eleshaky.
Software
sensitivity
Aerodynamic
anal-
Journal
of
design optimiza-
tion using sensitivity analysis and computational ics. AIAA Journal, 30(3):718-725, 1992.
shape
design
was first
[41 G.W. Burgreen and O. Baysal. Aerodynamic tion using preconditioned conjugate gradient paper
a finite
volume
formulation
to the degree configurations
and general
finite vol-
In the last year the technique
has been
to perform
that aerodynamic
design
with very rapid turnaround
can be accomplished
method. For demonstration ration with closely coupled
in a routine
July
fluid dynam-
shape optimizamethods. AIAA
1993.
fashion
and O. Baysal.
optimization
Three-dimensional
of wings
aerody-
using sensitivity
AIAA paper 94-0094, 32nd Aerospace Exhibit, Reno, Nevada, January 1994.
Sciences
has
of complete
I7] M. E. Eleshaky
is possible.
position
scheme
and O. Baysal.
analysis. In Proceedings fluid dynamics conference,
with the current
purposes a typical supersonic configunacelles and diverters was used as a test
[81 M. E. Eleshaky nacelle
case. Nacelles of the axisymmetric type were employed and situated such that they almost abutted the wing lower surface. The diverters
Preconditioned
for three-dimensional
and
Improving procedures.
domain
aerodynamic
decom-
sensitivity
of of the 12th AIAA computational pages 1055-1056, July 1993.
and O. Baysal.
Shape
optimization
near a flat plate wing using multiblock
ysis. AIAA paper 94-0160, 32nd Aerospace and Exhibit, Reno, Nevada, January 1994.
10
analysis.
Meeting
[61 G. W. Burgreen, O. Baysal, and M. E. Eleshaky. the efficiency of aerodynamic shape optimization AIAA paper 92-4697, September 1992.
implementa-
here, the technology shape
namic shape
the aerodynamic
[91. With the parallel
flow solver presented
93-3322,
[5] G. W. Burgreen
In this paper we have shown how the complicated design of supersonic aircraft configurations including airframe/nacelle interaction effects
will be
and flows governed
It has been demonstrated
mapping
participants
configurations
tion of the multiblock aircraft
[40].
by some industry
advanced
the techniques
meshes
ysis methods for the compressible Euler equations. Fluids Engineering, 113(4):681-688, 1991.
to perform calculations with arbitrary numerically generated grids [39, 25]. Further, results have been presented for three-dimensional
of future
has benefited
under
[1] O. Baysal
so as to
region with high Cps, thus pressure drag was reduced
[20, 39, 25, 23].
that it can be successfully
design
efforts,
13] O. Baysal and M. E. Eleshaky. Airfoil shape optimization using sensitivity analysis on viscous flow equations. Journal of Fluids Engineering, 115:75-84, 1993.
by the Euler equations
adopted
In future
References
that the
by the fourth author [19], the method has been verified by implementation for both potential flow and flows modeled
calculations
camber-only
and design cases were performed in less than one day, demonstrat-
siderable thanks also goes to David for his help in assembling the text.
CONCLUSIONS In the period
using the new
A preliminary
the cruise point performance for subject to wing thickness con-
both unstructured equations.
number N00014-92-J-1976, Research branch of NASA
by 3.5% or from CD = 0.00882 to CD = 0.00851. It is confidently believed that much more could be gained through the investment of more computer resources.
proposed numerical
were fully modeled method.
in Figure 8 are more
shocks impinge
a forward facing The final overall
surface
of the cross-sectional
shown
the nacelles
cusp the region where the nacelle surface. reducing
AFOSR
for both the upper and have been made to the
indicative of the design differences. The strong oblique shock ran near the leading edge has been largely softened. Meanwhile, lower surface
resources.
This research
run. Figure 7 shows the final solu-
field. The comparisons
lower
ACKNOWLEDGMENTS
result. It nevertheless demonstrates the power of a parallel implementation of the adjoint-based design method. The total wall clock
cuts of both the airfoils
and body
and design
extended to address by the Navier-Stokes
and time, only a two-design-
tion at the end of the second design iteration lower surfaces. It is clear that large differences
and the wing
ing that complete configuration designs may be achieved with rapid turn-around even with the most conservative estimates of available
cycle run of the optimization algorithm was achieved using 16 processors on an IBM SP2. This must be thought of as a very preliminary
time was 1.75 hours for the entire
wing
analysis
straints and fixed lift. All analysis on parallel architecture machines
The rest of the
complex mesh system in and around the nacelles shown was perturbed automatically as the design proceeded.
the nacelles
design was carried out to optimize the complete aircraft configuration
upper and lower surface for the baseline configuration at the cruise point. The nacelles and diverters were forced to follow the variaincluded
between
were constructed to eliminate the possibility of boundary layer entrainment into the nacelle inlets. These geometry entities as well
of a 3-D
sensitivity Sciences
anal-
Meeting
19]J. Gallman, Business tion.
J. Reuther, jet wing
AIAA paper
and Exhibit,
N. Pfeiffer, W. Forrest'
design
using
96-0554,
and D. Bemstorf.
aerodynamic
shape
34th Aerospace
Reno, Nevada,
January
[231
optimiza-
Sciences
A. Jameson. Theory, 1995.
Meeting
P. Gill,
W. Murray,
two revised mization. Numerical
and R. Pitfield.
quasi-Newton
The implementation
algorithms
R. M.
Hicks
and
optimization.
P. A. Henne.
[131 R. M. Hicks,
opti-
Division
of
January
J. Huan
and V. Modi.
1978.
and G. N. Vanderplaats.
of minimum
drag
and J.J. Alonso.
[261
An
bodies
AF/NAS
A. Jameson,
Bussoletti,
R.G. Melvin,
M.B. Bieterman,
D.P. Young,
F.T. Johnson,
and C.L. Hilmes.
Practical
and optimization in computational fluid dynamics. per 93-3111, AIAA 24th Fluid Dynamics Conference, Florida, July 1993.
127] R. Kennelly.
[17]
J.E. design
A. C. Taylor
Solution
A. Jameson.
Multigrid
algorithms
September
A. Jameson.
Aerodynamic
of Scientific
Computing,
1201 A. Jameson.
Automatic
design via control 3:233-260,
Mathe-
H. E. Jones.
design
Langley
[211 A. Jameson.
5-17,
February
Optimum
trol. In AGARD-VKI
airfoils
in Aerodynamics. 1994. A. Jameson. control
theory.
yon Karman
Optimum
Institute
aerodynamic
AIAA paper
tional Fluid Dynamics
design
to reduce
95-1729,
Conference,
via boundary
Optimum
Design
AIAA
San Diego,
using CFD
strategy.
fluid
HI, and P. A. Newman.
Aero-
using a 3-D supersonic of sensitivity derivatives.
Hampton,
Panama
NASA
consis-
TM 104207,
VA, February
S. Ta'asan, and M. D. Salas.
City,
G. W. Hou, and
for calculating
derivatives.
Euler A/AA
Symposium
and Optimization,
strategy
In
dynam&s
1992.
Airfoil optimization
1994.
and O. Baysal.
Design
optimization
of single-
and two-element airfoils on multiblock 4273, 5th AIAA/USAF/NASA/ISSMO
grids. AIAA paper 94Symposium on Multi-
disciplinary September
Panama
Analysis 1994.
and Optimization,
City, Florida,
133] J. Lewis and R. Agarwal. Airfoil design via control theory using the full-potential and Euler equations. Technical report,
con-
Methods
for Fluid Dynamics,
design
Center,
for Fluid Dynamics,
Journal
1990.
Series,
Research
1321 J.M. Lacasse
of transonic
iterative
by the one-shot method. In AGARD- VKI Lecture Series, Optimum Design Methods in Aerodynamics. von Karman Institute
1988.
aerodynamic Lecture
CFD sensitivity
1311 G. Kuruvila,
Methods,
the shock induced pressure drag. In Proceedings of the 31st Israel Annual Conference on Aviation and Aeronautics, Tel Aviv, pages
An incremental
tent discrete
E A. Newman,
for three-dimensional
5th AIAA/USAF/NASA/ISSMO Analysis 1994.
January
Applied AerodyJuly 1983.
1301 V. M. Korivi, A. C. Taylor HI, E A. Newman,
flow calcu-
theory.
94-4270,
on Multidisciplinary Florida, September
for two dimen-
for compressible
airfoil design-by-
of of the 12th AIAA computational pages 1053-1054, July 1993.
A. Jameson.
of the Euler equations
for transonic
Euler code using incremental
paper
91-3083,
solutions
with Runge-
81-1259,
supersonic
A. C. Taylor
de5th
AIAA paper
Proceedings conference,
viscous
AIAA paper
Numerical
derivatives
and external 1991.
flows.
optimiza-
FL, Septem-
methods
ffl, G. W. Hou,
Sensitivity
dynamic optimization studies code with efficient calculation
ings of the 2rid European Conference on Multigrid Cologne, 1985, Springer-Verlag, 1986.
122]
method
A.C. Taylor III, G. W. Hou, and V. M. Korivi. Sensitivity analysis, approximate analysis, and design optimization for internal
lations. In W. Hackbusch and U. Trottenberg, editors, Lecture Notes in Mathematics, Vol. 1228, pages 166-201. Proceed-
[19]
schemes.
Improved
and H. E. Jones.
AIAA paOrlando,
sional transonic flow by a multigrid method. Applied matics and Computations, 13:327-356, 1983. [18]
and E. Turkel.
optimization. AIAA paper 83-1864, AIAA namics Conference, Danvers, Massachusetts,
1291 V. M. Korivi, [161
Wiley,
on Muitidisciplinary
City Beach,
by finite volume
Kutta time stepping 1981.
in
Panama
W. Schmidt,
[281 V. M. Korivi, W.P. Huffman,
aerodynamic
A/IS SMO Symposium
of the Euler equations
incompressible laminar flow. Technical report, The Forum on CFD for Design and Optimization, (IMECE 95), San Francisco, California, November 1995. [15]
Automatic
Analysis and Optimization, ber 1994.
by numerical
by numerical optimization. NASA Center, Moffett Field, California,
Design
Control
1995.
1251 A. Jameson and J. Reuther. Control theory based airfoil sign using the Euler equations. AIAA paper 94-4272,
July 1974. [14]
Using
Review
1996.
AIAA/US design
15:407-412,
E. M. Murman,
assessment of airfoil design TM X-3092, Ames Research
Design
tion on distributed memory architectures. AIAA paper 96-0409, 34th Aerospace Sciences Meeting and Exhibit, Reno, Nevada,
PracticalOptimization.
Wing
JournalofAircrafi,
of
for unconstrained
NAC 11, National Physical Laboratory, Analysis and Computing, 1972.
[11 ] P.E. Gill, W. Murray, and M.H. Wright. Academic Press, 1981. [12]
Aerodynamic Fluid Dynam&s
1996. [241 A. Jameson
[10]
Optimum
Computational
The Forum on CFD for Design and Optimization, San Francisco, California, November 1995.
and
[341
12th Computa-
J. L. Lions.
Optimal
Control
tial Differential Equations. Translated by S.K. Mitter.
CA, June 1995.
11
of Systems
Springer-Verlag,
(IMECE
Governed
95),
by Par-
New York, 1971.
[351P.A.Newman, V. M. Korivi.
G. W. Hou,
for use in large-scale incorporating ley Research [361
137]
A. C. Taylor
on computational
gradient-based,
Ill, and
multidisciplinary
Optimal Shape Design New York, 1984.
O. Pironneau.
Optimal
AGARD-VKI
Lecture
Aerodynamics.
shape Series,
yon Karman
S. Cliff,
for
design
design
R. Hicks,
Elliptic
Design
[49]
In
Methods
and C.P. van Dam.
in
1994.
Practical
de-
using the Euler
[391 J. Reuther and A. Jameson. Control theory based airfoil design for potential flow and a finite volume discretization. AIAA paper 94-0499, 32rid Aerospace Reno, Nevada, January 1994. [40]
J. Reuther
and A. Jameson.
Sciences
Aerodynamic
of wing and wing-body configurations AIAA paper 95-0123, 33rd Aerospace Exhibit, [411
Reno, Nevada,
January
J. Reuther
and A. Jameson.
for control
theory
6th International
based
Meeting
1995.
A comparison
ics, Lake Tahoe, Nevada,
optimization
using control theory. Sciences Meeting and
of design
airfoil optimization.
Symposium
and Exhibit,
shape
on Computational
September
variables
Technical
report,
Fluid Dynam-
1995.
I421 J. Reuther and A. Jameson. Supersonic wing and wing-body shape optimization using an adjoint formulation. Technical report, The Forum on CFD for Design and Optimization, 95), San Francisco, California, November 1995. [431 J. Reuther, A. Jameson, ders. Aerodynamic figurations
via an adjoint
34th Aerospace January 1996. [44]
J. Reuther,
J. Farmer, L. Martinelli,
shape optimization Sciences
formulation. Meeting
airfoils with reduced 29:297-298, 1992.
[45]
J. J. Reuther. Aerodynamic theory. Ph.D. Dissertation, Davis, CA, June 1996.
[46]
J.P. Steinbrenner, GEN report,
Flight
velopment 1990. [471
S. Ta'asan,
block
Dynamics
Center,
Laboratory,
G. Kuruvila,
Aerospace Sciences uary 1992.
and C.L. Fouts.
grid generation
Wright-Patterson
sign and optimization
and transonic
pitching
moments.
shape optimization using conlrol University of California, Davis,
J.R. Chawner,
3D multiple
96-0094,
Reno, Nevada,
C.P. van Dam, and R. Hicks. Subsonic
low-Reynolds-number Journal of Aircraft,
aircraft con-
AIAA paper
and Exhibit,
(IMECE
and D. Saun-
of complex
Wright
Meeting
Technical
Research
and De-
Air Force Base, Ohio, July
and M. D. Salas.
in one shot.
The GRID-
system.
Aerodynamic
AIAA paper
and Exhibit,
Reno,
Z.U.A
Grid Generation, ence Publishing
Foundations and Applications. Elsevier Company, New York, NY, 1985.
Warsi,
and C.W.
D.P. Young,
W.P. Huffman,
R.G. Melvin,
C.L. Hilmes, and F.T. Johnson. vergence in design optimization.
Systems.
for Fluid Dynamics,
sign optimization of wing/body configurations equations. AIAA paper 92-2633, 1992.
J.F. Thompson,
Mastin.
Numerical Sci-
Lang-
for aerodynamics.
Optimum
Institute
[48]
methodologies
advanced CFD codes. NASA TM 104206, Center, Hampton, VA, February 1992.
O. Pironneau. Springer-Verlag,
[381 J. Reuther,
H. E. Jones,
Observations
de-
92-0025,
30th
Nevada,
Jan-
12
AIAA/US
AF/NASA/ISSMO
Analysis 1994.
and Optimization,
M.B.
Bieterman,
Inexactness and global conAIAA paper 94-4286, 5th
Symposium Panama
City,
on Multidisciplinary Florida,
September
Block
II
Center
! Block I
I
, , i
m , ,
,i ,r: i
I
i , i J.,.J...J i
I
_:
.-"
i
_ :i ,
i
Including
/
]
"1" I "1"
"1" I -I"
I
I
"1 I -I I
Solid
Block I
-
Center
i_"i'.'"'"_
/ Block
"
i --, i ] ] ] .......i .......i .......i,,Z,,.,i
III
Double
Boundary
Halo
I
Center
Figure 3:4 Block interface using Each block's double halo of cells
a double contains
adjacent
halo of cells around each block. information from internal cells in
blocks.
13
!
tin
0
I
II
||
("4
,r"4
o_
I
l I
_
_.Jf'_ _
It
16
.tmt
I
Q_
QI'3
,-,
_
it
N o_
o_
0
17
_
_
inilill_tan
8a: span station
_
_
)1 = 0.116
inkial
8b: span station
_okNiolt
_
_
7/= 0.374
inRial
Iokltiott
it
8c: span
station
z = 0.632
8d: span
Figure 8: SYN87-MB Fixed Lift Drag Minimization. 180-Block Mesh, 1,500 K mesh cells, M = 2.2 90 Camber ...... --,
Hicks-Henne
variables.
, Initial Pressure
Coefficient
Final
Pressure
18
Coefficient.
station
z = 0.871