Aerodynamic Shape Optimization of Supersonic Aircraft ... - CiteSeerX

3 downloads 0 Views 1MB Size Report
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,



+ _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



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