NOTICE

• J .F;_"

NASA CASE NO. LAR 14815-1-CU

PRINT FIG. 1

The invention disclosed in this document resulted from research in aeronautical and

space activities performed under programs of the National Aeronautics and Space

Administration. The invention is owned by NASA and is, therefore, available for licensing

in accordance with the NASA Patent Licensing Regulation (14 Code of Federal

Regulations 1245.2).

To encourage commercial utilization of NASA-owned inventions, it is NASA policy to grant

licenses to commercial concerns. Although NASA encourages nonexclusive licensing to

promote competition and achieve the widest possible utilization, NASA will provide the

necessary incentive to the licensee to achieve early practical application of the invention.

Address inquiries and all applications for license for this invention to NASA Patent

Counsel, Langley Research Center, Code 143, Hampton, Virginia 23665-5225.

Serial No.: 07/820,432

01-03-92 LaRC

(NASA-Case-LAR-I48IS-l-CU) AEROOYNAMIC

DESIGN OPTIMIZAT[ON USING SENSITIVITY

A_ALYSIS AN9 COMPUTATIONAL FLUID DYNAMICS

Patent Application (NASA) 24 p

G3/36

N92-29830

Unclas

0111101
NASA Case No. LAR 14815-1-CU 

AWARDS ABSTRACT 

AERODYNAMIC DESIGN OPTIMIZATION USING SENSITIVITY ANALYSIS 

AND COMPUTATIONAL FLUID DYNAMICS 

An efficient aerodynamic shape optimization method based on a 

computational fluid dynamics/sensitivity analysis algorithm has 

been developed which determines automatically the geometrical 

definition of an optimal surface starting from any initial 

arbitrary geometry. This method is not limited to any number of 

design variables or to any class of surfaces for the shape 

definition. 

The novelty of this invention is found in determining the 

optimal aerodynamic surface for a given set of parameters and 

constraints using a minimum number of computational fluid dynamics 

analyses. 

Inventors: Oktay Baysal Mohamad E. Eleshaky 

Address : 1346 Brunswick Avenue 18-D Elizabeth Road 

Norfolk VA 23508-1345 Hamp 

ss: - ss: 

Employer: Old Domion University Research Foundation 

Evaluator: David S. Miller 

Serial No.: 07/820, 432 

Filed: January 3, 1992 
(_j " .,

LAR 14815-1-CU -1 - PATENT APPLICATION

AERODYNAMIC DESIGN OPTIMIZATION USING SENSITIVITY ANALYSIS

AND COMPUTATIONAL FLUID DYNAMICS

5

Origin of the Invention

The Invention described herein was made in the performance of work

under NASA Grant No. NAG-l-1188. In accordance with 35 U.S.C. 202, the

grantee elected not to retain title.

10

Background of the Invention

1.. Technical Field of the Invention

15

The present invention relates generally to designing aerodynamic

shapes, and more particularly to the optimization of aerodynamic shapes using

sensitivity analysis and computational fluid dynamics.

2_=. Discussion of the Related Art

20

25

30

Aerodynamic shape optimization has recently become a subject of

considerable interest. This field of study involves the ability to determine the

geometry of an aerodynamic configuration that achieves specified objectives

subject to certain constraints and satisfies the governing equations of its

flowfield. For instance, in attaining the geometry of an optimum airfoil section

two basic categories of design methods exist.

The first category includes the inverse design methods, in which one

specifies the surface pressure distribution and the method calculates the

corresponding airfoil geometry. Examples of this approach include the

hodograph methods which solve the linear potential equations in the hodograph

plane, and a method which treats the two-dimensional flow as a set of

streamtubes coupled through the position of, and pressure at, the stream
5

10

15

20

25

LAR 14815-1-CU -2- PATENT APPLICATION

Interfaces. The latter method solves the Euler equations which guarantees a

proper treatment of shocks. Some of these methods have recently included the

viscous effects through a coupled integral boundary layer analysis.

The second category includes the numerical optimization methods, which

couple a fiowfield analysis algorithm with an optimization algorithm. The

aerodynamic quantities such as lift, drag, pitching moment, and pressure

distribution are computed by the flowfield analysis algorithm for a certain

configuration and are used in defining an objective function to be minimized or

maximized. This objective function must relate changes in geometry to

improvements in the aerodynamic quality of the design.

These methods require a significant level of experience and skill to

define the geometrical aerodynamic shapes.

It is accordingly an object of the present invention to design an

aerodynamic shape which maximizes performance within specified aerodynamic

and geometrical constraints.

It is another object of the present invention to determine the unique

aerodynamic shape which maximizes performance within specified aerodynamic

and geometrical constraints regardless of the initial shape chosen to optimize.

It is another object of the present invention to achieve the foregoing

objects with a process which does not require a complete computational fluid

dynamics analysis of the shape as part of every iteration.

It is yet another object of the present invention to accomplish the

foregoing objects in a simple manner.

Additional objects and advantages of the present invention are apparent

from the drawings and specification which follow.

3O

Summary of the Invention

According to the present invention, the foregoing and additional objects

are obtained by choosing a function to optimize, determining a set of physical

constraints, a set of constant design parameters and a set of variables which
5

10

15

20

25

LAR 14815-1-CU -3- PATENT APPLICATION

may be modified to optimize the function. An initial shape must also be

chosen, although the result is not dependent on this initial shape as this

process determines the unique optimal shape. The constraints are analyzed

to determine what the feasible space is, i.e. what ranges must each variable is

limited to and an Initial value for each variable is chosen within this feasible

space. A complete computational fluid dynamics analysis is performed for the

initial shape, the initial set of constant design parameters and the initial set of

variables. The sensitivity coefficients for the initial shape and variables are

calculated and the function is evaluated for the initial set of variables, ensuring

that none of the constraints are violated. Based on the evaluation of the

function, it is determined whether one or more than one variable will be

changed. If one variable is changed, it is changed less than about 5% and the

shape is redefined based on this change. An approximate flow analysis for the

redefined shape is performed based on the sensitivity coefficients previously

calculated. The function is evaluated for the new variables, ensuring that none

of the constraints are violated and the redefined shape is evaluated to

determine if the shape is optimal. If more than one variable is changed, the

shape is redefined based on these changes and a complete computational fluid

dynamics analysis is performed for the redefined shape. The function is

evaluated for the new variables, ensuring that none of the constraints are

violated and the redefined shape is evaluated to determine if the shape is

optimal. These steps are repeated until the shape is optimal. Once the optimal

shape is obtained, the shape may be evaluated to determine the sensitivity of

the shape with respect to one or more of the design parameters.

Brief Description of the Drawings

Fig. 1 is a flow chart showing the overall design optimization process;

Fig. 2 is a flow chart showing the details of the direct method of the

sensitivity analysis process; and
LAR 14815-1 -CU -4- PATENT APPLICATION

Fig. 3 is a flow chart showing the details of the adjoint variable method

of the sensitivity analysis process.

5

Detailed Description of the Invention

The overall process of aerodynamic design optimization consists of the

computational fluid dynamics analysis, the sensitivity analysis, the approximate

flow analysis and the optimization method. A flow chart of the process is

presented in Fig. 1. The following example illustrates this process.

10

15

Example

It is desired to determine the angles of the nozzle ramp, o_,and of the

cowl, 13, that yield a maximum axial thrust force coefficient, F, subject to

constraints, Gj. The angles 0¢and 13are the design variables, X_ and therefore,

the number of design variables, NDV, is equal to two.

Mathematically, it is required to get

maxF(Q(X@,X

(1)

2O

Gj(O(X_,X_ s o;

Xo,.., s X o < Xo,_

j= 1,NCON (2)

(3)

where F is the objective function, NCON is the number of constraints, and Q

is the vector of the conserved variables of the fluid flow. XD,o_r and Xo

upper
5

LAR 14815-1 -CU -5- PATENT APPLICATION

are the lower and the upper bounds of the design variables. A judicious

choice of the upper and lower bounds for the design variables accelerates

convergence of the optimum solution. A fixed nozzle length, a fixed cowl

thickness, and a fixed nozzle inlet height are imposed as the geometric

constraints.

The axial component of the thrust force due to nozzle wall shape,

F,=_, is obtained by integrating numerically the pressure over the ramp and

cowl surfaces.

10 (4)

This force is normalized by the force associated with the inflow given by

F_ n, - [; P_(I+'yM2_dY

15

(5)

In the case of an inflow parallel to the cowl with a constant Mach number,

this force is centered at the mid-point of the line segment k--c,and its value is

20

By definition, the axial thrust force coefficient is given by

(6)
LAR 14815-1 -CU -6- PATENT APPLICATION

F m

F_dm/

F/n/Ibw

(7)

5

This axial thrust force coefficient is subject to the following three constraints

(i.e., NCON = 3):

1. The static pressure at the ramp tip, P_,is forced to reach a

percentage C1, of the freestream static pressure, P., such that a maximum

expansion is reached without any back flow, i.e.

G,(X_ - 1 PI _ o

c,P.

10

15

(8)

2. The static pressure at the cowl tip, Pn, should be within

specified limits, (C2 and C3), of the freestream static pressure, P., such that

a maximum expansion is reached without any back flow on either of the

cowl surfaces, i.e.

G2(X"_ - 1 pn .<0

C,P.

(9)

1sO

(10)
LAR 14815-1 -CU -7- PATENT APPLICATION

The derivatives of the objective function, F, and constraints Gp with

respect to the design variables, Xo, are given by

5

v'="

(11)

dX'o aX'. + kaO) J - 1, NCON

(12)

10

15

20

Two different quasi-analytical approaches are available for obtaining

these derivatives, the direct method, shown in figure 2, and the adjoint

variable method, shown in figure 3. The choice of the particular sensitivity

analysis approach depends on the number of design variables (NDV) and

adjoint vectors (NCON+I). If the number of design variables is less than the

number of adjoint vectors, the direct method is more efficient than the adjoint

variable method. In the present example, the number of design variables

(NDV = 2) is less than the number of adjoint vectors (NCON+I = 4), the

direct method is selected to determine the sensitivity coefficients.

The governing equations for two-dimensional, steady, compressible,

inviscid flow of an ideal gas with constant specific heat ratios written in the

residual vector form are
LAR 14815-1-CU -8- PATENT APPLICATION

5

a_/÷ o

a_ an

(13)

where f and g are the flux Jacobians in generalized coordinates (_, TI).

The quasi-analytical direct method begins with the differentiation of

Equation (13) with respect to the design variables to yield the sensitivity

equation:

(14)

10 Then Equation (14) is solved for o_/onXo. Equation (14) needs to be

I

solved once for each design variable Xo, however, the coefficient matrix

I !

(o_R/onQ)needs to be factorized once and for all. The remaining partial

derivatives in Equations (11) and (12) are evaluated analytically using

Equations (7-10). The final step is determining the values of VF and VG v

15 In order to reduce the number of CFD analyses during the optimization

process, an approximate flow analysis is performed, which is based on a

Taylor-series expansion of the vector of conserved variables Q(X o + _Xo) about

Q(X;) as follows:
LAR 14815-1-CU -9- PATENT APPLICATION

a(X_+ AXo)- (_'(X" _) + /8-_O) A_VO+...

x'o-,,'r_

(15)

Substituting Equation (14)into Equation (15) results in,

5

10

15

20

an-(a(x;),x;)] a'=r(cqx;)'xt;) Ao- - AXo

(16)

where AQ = Q(X; + AXo)- Q(Xo).

Equation (16) gives the changes in Q due to the changes in AX o. In

other words, given the flow field solution, Q(Xo), associated with a

configuration, XD, the flow field solution Q(X* + AXo), associated with the

configuration, (X_ + AXo), is obtained via Equation (16).

Although, the approximate analysis is less accurate than the actual

analysis, it is less costly in terms of computer time, especially, when the

number of flow field governing equations is large. However, the approximate

analysis is of acceptable accuracy up to +5% changes in the design

variables.

Equation (13) is solved by a first-order accurate, implicit, upwind finite

volume scheme. An upwind discretized form of Equation (13) at an interior

cell (i,j)is as follows:
LAR 14815-1-CU -10- PATENT APPLICATION

_1_1"_((_-l+112_Mi+112) - _+((_-/-1/2_M1-1/2,/_

+ _t(_+/,lrj_Ml+1/2_- _-(_+l_112jbMl_l12_

÷ t(O-_,:,M,j,,=) - _(O-u_,:,Mu_,_)

+ _-(O'_/+I/2,MIj.,lt,_- _-(0"_1_112,M1__112 )

(17)

5

where f+, f', g+, and g" are the operator-split inviscid fluxes. M represents

the projected surface areas for these fluxes and it is associated with the

--+

coordinate transformation metrics. The Q i+l_,j and Q_÷I_.jfor a first-order

accurate scheme are given by,

(18)

10 Consequently, the left-hand side of Equation (14) is evaluated by

differentiating the upwind discretized form given by Equation (17) for an

Interior cell (i,j) with respect to (Q:1:_.1_,l_112) as follows,

[{.._ } {__X_t]lj =, ¢_jl_l_/ ¢_(_-1÷1/2,_ + ¢_j_i;l_/ a(_-/_l/p. _

(TN_-h-1/P.J c_X D 0_0-/- 1/2,/ ¢3X D

* aO÷_÷lr_jaXo aO÷_-lr_j eX.

a,qq aO-,_,_= a,_v aO-,,,__

+ 80- .lr_ aXD + aO_l,/_lr 2 (3XD

aRq aQ÷q.iR + aR v aQ÷v_I=

* C_(_÷lj,112 aXD C'/0÷1,/..1/2 aX D

(19)
LAR 14815-1 -CU -11 - PATENT APPLICATION

Denoting the flux Jacobians by A and B, Equation (19) can be rewritten as,

- A+ u a_ °

+

8)( ° - A-i/ 8)( °

- aO-,t_lr

+ B,uSQ-I,/*ll 2 _ B÷i,l_l

8)(o 8Xe

* B-U+1 8X ° - U 8X °

(20)

5

10

The values of the conserved variables at the cell-interface locations (Q*P.1_,_1_2)

can be viewed as functions of the values of the conserved variables of the

neighboring cells evaluated at the centers Therefore, the terms [o_l:_.,_j/o_Xj

and [o_+j,l_o_Xo] in Eq. (20) can be expressed in terms of [c3Q_,t'o%XD],

[o_Q,Jo_,Xj,and [o_Q,j,,/c3Xo].Hence, Eq (20)is rewritten as

- A+I_1 + (A÷¢B'¢A-¢B-,._aQ-__ I

q 8Xo

+ A-,.I_+B+...SQ-__-_+B-...SQ___ t _

" aX o ,.r. aX o ,a+, aXe

(21)

15

n D

The coefficients A through E in Equation (21) are 4x4 blocks which are

functions of the flux Jacobians A" and B+. These flux Jacobians A* and B* are

usually available for an interior cell from the flowfield solution when forcing the

residual to zero. Hence, only a few computations are needed to assemble the
5

10

15

LAR 14815-1-CU

w

Jacobian matrix [_R/o-_)] in Eq (14).

residual expression (Equation 17)

boundary conditions.

-12- PATENT APPLICATION

However, it is necessary to revise the

at the boundary points to include the

Once the numerical values of all the elements in the (n x n) coefficient

matrix are obtained, the coefficient matrix can be assembled. It has a block

banded structure with five nonzero computational domain with I cells in the P_

direction and J cells in the 11direction, the matrix dimension, n, is (41J). The

first order upwind discretization of the governing equations yields a coefficient

matrix with a subdiagonal bandwidth of (4J + 3) or (41 + 3), depending on the

way in which the unknowns are ordered. The diagonal storage of band

matrices and the unknown ordering resulting in the minimum bandwidth are

used in the present example to reduce the computer memory requirement.

The right-hand side of Equation (14) is evaluated by differentiating

Equation (17) for an interior cell (i,j) with respect to the design variables, Xo, as

follows,

2O

-I-

)

a(M,_,,)

8¢ (O-,,/-,p, Mq-, R)

+

6(M,.,/'z,_

+

_M,-,rz)

+ a¢(('_+q.'12'Mlp'l'z)

a(Mu_, 

(22)

B

As it is seen from Equation (22), [c3R/_Xo] depends on the derivatives of

the projected surface areas (M) with respect to the design variables. If an

analytical expression for M = M(XD) exists, then this differentiation is
5

10

15

2O

25

LAR 14815-1-CU -13- PATENT APPLICATION

straightforward. Otherwise, a finite difference approximation for o_M/o_Xowith a

small step size AXD can be used.

Finally, the solution of the system of linear algebraic equations (14) is

achieved by using the Gauss elimination method on the elements within the

bandwidth only.

The computational fluid dynamics analyses are performed using the

general purpose, finite volume Euler/Navier-Stokes CFD code "VUMXZ3". This

code has been applied to analyze a variety of complex internal and external

flows. It produces consistent and repeatable flow simulations in the sense that

small perturbations to design variables are accurately reflected in the flowfield

solution.

A key part of the present design procedure is the sensitivity analysis

where the derivatives of the constraints, the objective function and the

conserved flow variables, with respect to the design variables are computed.

The derivatives quantify the effects of each design variable on the design and

thereby identify the most important design changes to make enroute to the

optimum design. Two different quasi-analytical approaches are available for

obtaining these derivatives, the direct method, shown in figure 2 and the adjoint

variable method, shown in figure 3. The choice of the particular sensitivity

analysis approach depends on the number of design variables and adjoint

vectors. If the number of design variables is greater than the number of adjoint

vectors, the adjoint variable method is more efficient than the direct method.

The present nonlinear, constrained optimization problem is solved using

the modified feasible directions method developed by Vanderplaats ("An

Efficient Feasible Direction Algorithm for Design Synthesis," AIAA J., Vol. 22,

No. 11, October, 1984, pp. 1633-1640). Given a set of initial conditions for

design variables, the method first determines if the initial values are in the

feasible design space, i.e., the space when none of the design constraints are

violated. If not, then the design variables must be changed until the feasible
LAR 14815-1-CU -14- PATENT APPLICATION

design space is located. At this point, the method can search to find the

optimum design within the feasible design space.

Since the optimization process requires many evaluations of the

objective function and constraints before an optimum design is obtained, the

5 process can be very expensive if a CFD analysis were performed for each

evaluation. However, the optimization process primarily uses analysis results

to move in the direction of the optimum design. Hence an analysis needs to

be made only occasionally during the design process and always at the end to

check the final design. In the present process, approximate flow analyses

10 using Equation 16 are performed through the one-dimensional search of the

optimization process, whereas the CFD analysis is performed only when new

gradients of constraints and objective function are needed, i.e. when the design

changes substantially.

Finally the sensitivity information is passed to the optimizer along with

15 the current values of the design variables, constraints, and the objective

function. The optimizer uses this information to generate a new set of design

variables and the entire procedure is repeated until a converged design is

obtained. The design algorithm is deemed converged when all the constraints

are satisfied and the objective function has a value which does not change for

20 a specified number of optimization loops.

After the optimum design is obtained, it is desirable to determine the

sensitivity of the optimum design with respect to one or more design

parameters. Such information is useful to perform trade-off analyses. For

example, it may be wished to estimate what effect a specified increase in the

25 free stream Mach number has on the optimum thrust. Mathematically, this

requires the derivatives of the optimum values of the objective function and the

corresponding design variables with respect to the design parameters. In the

present example, the design parameters are the specific heat ratio, 7, Mach

numbers, total temperatures and total pressures of both the external flow and

30 the internal flow leading to the nozzle jet.
5

LAR 14815-1 -CU -15- PATENT APPLICATION

The first-order sensitivity derivative method is adapted for the present

invention. The vector P contains the design parameters, which are held fixed

during the optimization. Using the superscript "op" to denote the optimum

quantities, the dependence of F°p and G on Xo and P can be written as,

FOp. I_(Q(X D(t_),t_),X o(P),l _)

(23)

_,- G a(Q(X o(P),P),X D(P),P) - 0

(24)

10 where G. is a vector containing only the active constraints at the constrained

maximum. The total optimum sensitivity derivative of the objective function with

respect to a design parameter P is obtained using the chain rule of

differentiation

15 (25)

Any perturbation of the parameter P about its value at the initial optimum must

be such that the originally active constraints remain active,
LAR 14815-1-CU -16- PATENT APPLICATION

(26)

The Kuhn-Tucker condition satisfied at a local optimum is

5 (27)

where _ is a vector containing the Lagrangian multipliers. Combining Eqs. (26)

and (27) yields,

10

i

Equation (28) is put into Equation (25) in order to eliminate i_--_XI'

_P /

(28)

(29)
LAR 14815-1-CU -17- PATENT APPLICATION

The Lagrangian multipliers in Equation (29) are obtained from the following

relation evaluated at the optimum point

5

-

(30)

m

The derivatives of the vector of conserved variables Q with respect to the

design parameters are obtained as follows

10

t_(U(Xo(P),P),Xo(P),P)- 0

(31)

Differentiation of Equation (31) with respect to the design parameters yields,

(32)

15 Using Equation (14) with Equation (32) results in,

20

(33)

Solving Equation (33) similarly to Equation (14) for _ and substituting it into

(29) yields the sensitivities of the objective function to the design parameters.
5

LAR 14815-1 -CU -18- PATENT APPLICATION

Since the number of design parameters (equal to seven to the present

example) is greater than the number of the adjoint vectors (equal to one for the

present example), the adjoint variable method, is more economical. This can

be seen from the following substitution of Equation (33) into Equation (29).

dF °p . _aF°p . "_T _aG- m + Ifa l T

dP aP aP [_ aQ /

(34)

10

where j-1- [_Q]-I

equation is defined,

Then, an adjoint vector t, that satisfies the following

JT'_ " [ aF_ + (_-_" aGla _ )]_

(35)

Substitution of Equation (35) into Equation (34) gives

15

dr-'_

dP aP aP

(36)

The adjoint system of Equation (35) is independent of any differentiation

with respect to the design parameters. Also, both terms on the right hand side

of Equation (35) are available from the calculations of Eqs. (11) and (12). The
LAR 14815-1-CU -19- PATENT APPLICATION

partial derivatives, _FOp

o_p ._p and-_p are evaluated analytically. Therefore, the ,

sensitivity derivatives (Equation 36) are obtained after solving Equation (35),

evaluating Equation (30) and pertinent substitutions.

What is claimed is:
LAR 14815-1 -CU PATENT APPLICATION

AERODYNAMIC DESIGN OPTIMIZATION USING SENSITIVITY ANALYSIS

AND COMPUTATIONAL FLUID DYNAMICS

5

Abstract of the Disclosure

An efficient aerodynamic shape optimization method based on a

computational fluid dynamics/sensitivity analysis algorithm has been developed

which determines automatically the geometrical definition of an optimal surface

starting from any initial arbitrary geometry. This method is not limited to any

number of design variables or to any class of surfaces for the shape definition.
IOF3

AERODYNAMIC DESIGN OPTIMIZATION FLOWCHART

Problem Parameters Vor,ob,.s P J__Definition Sur'oce I I

Evaluate J Optimizer (ADS) I

F'Gj's "]' I

Approximate

Flow

Analysis

I New XD_

(=initial XD

in first pass)

I

I Definition Surface

INO

Flow Analysis

(CFD)

(VUMXZ3)

l

Computing Sensitivi t

Coefficients (Gradien[s)

VF, VG'j

using Sensitivity

Analysis Approach

YES

f

Computing the

Sensitivity of

the Optimum

Design w.r.t.P

FIG. 1
2OF3

I. DIRECT METHOD

I Flow Analysis

_____ i =k+l

J

_R eF aGi's aR

_' a_' _' _xo

Factorize (aR/_Q)

I _=' 1

¢

solve

I Compute, VF

_. Compute, VG j

NO i .

j=j., _ ,=_o. ]

¥ES

'YES

FIG. 2
30F3

13. ADJOINT VARIABLE METHOD +

I Flow Analysis j

__' a_' o_ ' axo

solve

(a_/aS )+. X_ = (aF/aS)

+

Compute, VF

j=l I

solve

FIG.3