I 

- 

NASA Contractor Report 181749 

ICASE REPORT NO. 88-64 

ICASE 

I 

I 

AERODYNAMIC DESIGN VIA CONTROL THEORY 

I 

Antony Jameson 

Contract No. NAS1-18107 

November 1988 

(bBSA-CR-lb1749) AEROCYYARIC UESIGD VIA 

CCYTBCL IRE06Y Biaal Beport IBASA) 42 p 

CSCL 01A 

G3/02 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 

National Aeronautics and 

Space Administration 

Hampton, Mrginia 23665 

~ROMWCticartw 

Y89-134C8 

Dnclas 

0179842 
AERODYNAMIC DESIGN VIA CONTROL THEORY 

Antony Jameson 

Department of Mechanical and Aerospace Engineering 

Princeton University 

Princeton, NJ 08540 

Abstract 

This paper addresses the question of how to modify in aerodynamic design to 

improve the performance. Representative examples are given to demonstrate the computational feasibility of using control theory for such a purpose. 

An introduction and historical survey of the subject is included. 

c 

This research was supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-18107 while the author was in residence at the Institute for Computer Applications in 

Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 

i 
1. Introduction and historical survey 

Computers have had a twofold impact on the science of aerodynamics. On the one hand 

numerical simulation may be used to gain new insights into the physics of complex flows. On the 

other hand computational methods can be used by engineers to predict the aerodynamic 

characteristics of alternative designs. Assuming that one has the ability to predict the 

performance, the question then arises of how to modify the design to improve the performance. 

This paper is addressed to that question. 

. 

Prior to 1960 computational methods were hardly used in aerodynamic analysis. The 

primary tool for the development of aerodynamic configurations was the wind tunnel. Shapes 

were tested and modifications selected in the light of pressure and force measurements together 

with flow visualization techniques. Computational methods are now quite widely accepted in the 

aircraft industry. This has been brought about by a combination of radical improvements in 

numerical algorithms and continuing advances in both speed and memory of computers. 

~ 

I 

1 

I 

I If a computational method is to be useful in the design process, it must be based on a 

mathematical model which provides an appropriate representation of the significant features of the 

flow, such as shock waves, vortices and boundary layers. The method must also be robust, not 

liable to fail when parameters are varied, and it must be able to treat useful configurations, 

I 

ultimately the complete aircraft. Finally reasonable accuracy should be attainable at reasonable 

cost. Much progress has been made in these directions [l-lo]. In many applications where the 

flow is unseparated, including designs for transonic flow with weak shock waves, useful predictions 

C 

I 

can be made quite inexpensively using the potential flow equation [14]. Methods are also 

available for solving the Euler equations for two- and three-dimensional configurations up to a 

complete aircraft [S-lo]. Viscous simulations are generally complicated by the need to allow for 

turbulence: while the Reynolds averaged equations can be solved by current methods, the results 

depend heavily on the choice of turbulence models. 
2 

Given the range of well proven methods now available, one can distinguish objectives for 

computational aerodynamics at several levels: 

1) Capability to predict the flow past an airplane or important components in different flight 

regimes such as take-off or cruise, and off design conditions such as flutter. 

2) Interactive design calculations to allow rapid improvement of the design. 

3) Automatic design optimization. 

Substantial progress has been made toward the first objective, and in relatively simple 

cases such as an airfoil or wing in inviscid flow, calculations can be performed fast enough that the 

second objective is within reach. The third objective has also been addressed for various special 

cases. In particular it has been recognized that the designer generally has an idea of the kind of 

pressure distribution that will lead to the desired performance. Thus it is useful to consider the 

problem of calculating the shape that will lead to a given pressure distribution. Such a shape does 

not necessarily exist, unless the pressure distribution satisfies certain constraints, and the problem 

must therefore be very carefully formulated: no shape exists, for example, for which stagnation 

pressure is attained over the entire surface. 

The problem of designing a two dimensional profile to attain a desired pressure distribution 

was first studied by Lighthill, who solved it for the case of incompressible flow by conformally 

mapping the profile to a unit circle [ll]. The speed over the profile is 

I where 4 is the potential for flow past a circle, and h is the modulus of the mapping function. The 

solution for 4 is known for incompressible flow. Let qd be the desired surface speed. Then the 

surface value of h can be obtained by setting q = qd in equation (l.l), and since the mapping 
3 

function is analytic, it is uniquely determined by the value of h on the boundary. A solution 

exists for a given speed q at infinity only if 

m 

and there are additional constraints on q if the profile is required to be closed. 

Lighthill's method was extended to compressible flow by McFadden [12]. Starting with a 

given shape, and a corresponding mapping function h('), the flow equations can be solved for the 

potential &'I, which now depends on h('). A new mapping function h(') is then determined by 

setting q = qd in equation (l.l), and the process is repeated. In the limiting case of zero Mach 

number the method reduces to Lighthill's method, and McFadden gives a proof that the iterations 

will converge for small Mach numbers. He also extends the method to treat transonic flow 

through the introduction of artificial viscosity to suppress the appearance of shock waves, which 

would cause the updated mapping function to be discontinuous. This difficulty can also be 

overcome by smoothmg the changes in the mapping function. Such an approach is used in a 

computer program written by the author for Grumman Aerospace. It allows the recovery of 

smooth profiles that generate flows containing shock waves, and it has been used to design 

improved blade sections for propellers [13]. A related method for three dimensional design was 

devised by Garabedian and McFadden [14]. In their scheme the steady potential flow solution is 

obtained by solving an artificial time dependent equation, and the surface is treated as a free 

boundary. This is shifted according to an auxiliary time dependent equation in such a way that 

the flow evolves toward the specified pressure distribution. 

Another way to formulate the problem of designing a profile for a given pressure 

distribution is to integrate the corresponding surface speed to obtain the surface potential. The 

potential flow equation is then solved with a Dirichlet boundary condition, and a shape correction 

is determined from the calculated normal velocity through the surface. This approach was first 
4 

tried ~y Tranen [15]. Volpe and Melnik have shown how to allow for the constraints that must be 

satisfied by the pressure distribution if a solution is to exist [IS]. The same idea has been used by 

Henne for three-dimensional design calculations [ 171. 

The hodograph transformation offers an alternative approach to the design of airfoils in 

transonic flows. Garabedian and Korn achieved a striking success in the design of airfoils to 

produce shock-free transonic flows by using the method of complex characteristics to solve the 

equations in the hodograph plane [MI. Another design procedure has been proposed by Giles, 

Drela and Thompkins [19], who write the two-dimensional Euler equations for inviscid flow in a 

streamline coordinate system, and use a Newton iteration. An option is then provided to treat the 

surface coordinates as unknowns, while the pressure is fixed. 

Finally, Hicks and Henne have explored the possibility of meeting desired design objectives 

by using constrained optimization [20]. The configuration is specified by a set of parameters, and 

any suitable computer program for flow analysis is used to evaluate the aerodynamic 

characteristics. The optimization method then selects values of these parameters that maximize 

some criterion of merit, such as the lift-to-drag ratio, subject to other constraints such as 

required wing thickness and volume. In principle this method allows the designer to specify any 

reasonable design objectives. The method becomes extremely expensive, however, as the number 

of parameters is increased, and its successful application in practice depends heavily on the choice 

of a parametric representation of the configuration. 

The purpose of this paper is to propose that there are benefits in regarding the design 

problem as a control problem in which the control is the shape of the boundary. A variety of 

alternative formulations of the design problem can then be treated systematically by using the 

mathematical theory for control of systems governed by partial differential equations [2l]. 

Suppose that the boundary is defined by a function f(x), where x is the position vector. As in the 

case of optimization theory applied to the design problem, the desired objective is specified by a 
5 

cost function I, which may, for example, measure the deviation from a desired surface pressure 

distribution, but could also represent other measures of performance such as lift and drag. The 

introduction of a cost function has the advantage that if the objective is unattainable, it is still 

possible to find a minimum of the cost function. Now a variation in the control M leads to a 

variation 61 in the cost. It is shown in the following sections that bl can be expressed to first order 

as an inner product of a gradient function g with 61: 

t 

61 = (g 1 M) 

Here g is independent of the particular variation 6f in the control, and can be determined by 

solving an adjoint equation. Now choose 

M = -xg 

where X is a sufficiently small positive number. Then 

61 = -X(g,g) < 0 

assuring a reduction in I. After making such a modification, the gradient can be recalculated and 

the process repeated to follow a path of steepest descent until a minimum is reached. In order to 

avoid violating constraints, such as a minimum acceptable wing thickness, the steps can be taken 

along the projection of the gradient into the allowable subspace of the control function. In this 

way one can devise design procedures which must necessarily converge at least to a local 

minimum, and which might be accelerated by the use of more sophisticated descent methods. 

While there is a possibility of more than one local minimum, the cost function can be chosen to 

reduce the likelihood of difficulties caused by such a contingency, and in any case the method will 

lead to an improvement over the initial design. The mathematical development resembles in 

many respects the method of calculating transonic potential flow proposed by Bristeau, 
6 

Pironneau, Glowinski, Periaux, Perrier and Poirier, who reformulated the solution of the flow 

equations as a least squares problem in control theory [4]. 

In order to illustrate the application of control theory to design problems in more detail, 

the following sections present design procedures for three examples. Section 2 discusses the design 

of two dimensional profiles for compressible potential flow when the profile is generated by 

conformal mapping. This leads to a generalization of the methods of Lighthill and McFadden. 

Section 3 discusses the same problem when the flow is governed by the inviscid Euler equations. 

Finally, Section 4 addresses the three dimensional design problem for a wing, assuming the flow to 

be governed by the inviscid Euler equations. The procedures which are presented require the 

solution of several partial differential equations at each step. The question of the most efficient 

discretization of these equations is deferred for future investigation. 
7 

2, Design for Dotential flow usina conformal mamine; 

Consider the case of two dimensional compressible inviscid flow. In the absence of shock 

waves an initially irrotational flow will remain irrotational, and we can assume that the velocity 

vector CJ is the gradient of a potential 4. In the presence of weak shock waves this remains a fairly 

good approximation. Let 5, T and S denote vorticity, temperature and entropy. Then according 

to Crocco’s Theorem, vorticity in steady flow is associated with entropy production through the 

relation 

gxL+ TV S = 0 

Thus, the introduction of a potential is consistent with the assumption of isentropic flow, and 

shock waves are modelled by isentropic jumps. Let p, p, c and M be the pressure, density, speed 

of sound and Mach number q/c. Then the potential flow equation is 

I 

V*pV@=O 

where the density is given by 

while 

Here M is the Mach number in the free stream, and the units have been chosen so that p and q al 
have the value unity in the far field. Equation (2.2) is a consequence of the energy equation in the 

form 

2 2 ,k +- 5 = constant 

D 

6 

z 

D 

Figure 1 

Suppose that the domain D exterior to the profile C in the z plane is conformally mapped 

onto the domain exterior to a unit circle in the u plane as sketched in Figure 1. Let R and B be 

polar coordinates in the u plane, and let r be the inverted radial coordinate 1/R. Also let h be the 

modulus of the derivative of the mapping function 

Now the potential flow equation becomes 

where the density is given by equation (2.1), and the circumferential and radial velocity 
9 

components are 

while 

22 q2=u +v 

The condition of flow tangency leads to the Neumann boundary condition 

v=:g=o on^ 

In the far field the potential is given by an asymptotic estimate, leading to a Dirichlet boundary 

condition at r = 0 [2]. 

Suppose that it is desired to achieve a specified velocity distribution qd on C. Introduce 

the cost function 

The design problem is now treated as a control problem where the control function is the mapping 

modulus h, which is to be chosen to minimize I subject to the constraints defined by the flow 

equations (2.2 - 2.7). 

A modification 6h to the mapping modulus will result in variations 64, SU, bv and bp to the 

potential, velocity components and density. The resulting variation in the cost will be 
10 

where on C q = u. Also 

while according to equations (2.2) and (2.6) 

Hence 

bp = -5 (u6u + v6v) 

C 

=P%T-%i(Ub4e+vrb4r) 26h 

C C 

It follows that bq5 satisfies 

a 26h a La4 = - 'i78 (PM 48 3) - r 5 (PM2r4r 4 

where 

Then if $ is any periodic differentiable function which vanishes in the far field 

% L64 dS = 5 pM2 V$ V$T 6h dS 

D r 

(2.10) 

(2.11) 

! where dS is the area element rdrd8, and the right hand side has been integrated by parts. 

I 
11 

Now we can augment equation (2.9) by subtracting the constraint (2.11). The auxiliary 

function $J then plays the role of a Lagrange multiplier. Substituting for bq and integrating the 

term 

by parts, we obtain 

Now suppose that $J satisfies the adjoint equation 

L$=O inD 

with the boundary condition 

Then integrating by parts 

and 

(2.12) 

(2.13) 

(2.14) 
12 

Here the first term represents the direct effect of the change in the metric, while the area integral 

represents a correction for the effect of compressibility. 

Equation (2.14) can be further simplified to represent bI purely as a boundary integral 

because the mapping function is fully determined by the value of its modulus on the boundary. 

Set 

where 

f = log lgl = log h 

and 

Then f satisfies Laplace’s equation 

Af=O inD 

and if there is no stretching in the far field, f + 0. Thus 

AM=O inD 

and M-, 0 in the far field. 

Introduce another auxiliary function P which satisfies 

2 AP=pM V#*V$ inD (2.15) 
13 

1 and 

P=O onC 

Then the area integral in equation (2.14) is I 

and finally 

J APbfdS = 1 HEdBD C 

&=JgbfdO 

C 

(2.16) 

(2.17) 

(2.18) 

This suggests setting 

6f= -xg 

so that if X is a sufficiently small positive number 

Arbitrary variations 6f cannot, however, be admitted. The condition that f + 0 in the far 

field, and also the requirement that the profile should be closed, imply constraints which must be 
14 

satisfied by f on the boundary C. Suppose that log[$] is expanded as a power series 

m 

log [%I = 1 > 

U n=O 

(2.19) 

where only negative powers are retained uecause otherwise $ would become unbounded for large 

u. The condition that f + 0 as u + m implies 

co = 0 

Also the change in z on integration around a circuit is 

Az=f.. dz da = 27ri c1 

so the profile will be closed only if 

c1 = 0. 

On C equation (2.19) reduces to 

m m 

fc + i pc = 1 (an cos n6 + bn sin n6) + i 1 (bn cos n8- an sin n 6) 

n=O n=O 

Thus an and bn are the Fourier coefficients of fc, and these constraints reduce to 

a. = 0, a, = 0, b, = 0 

In order to satisfy these constraints we can project g on to the admissible subspace for fc 
15 

where 

1 by setting 

= g - A. - A1 cos B-B1 sin 0 

A0 = & J 

C 

AI = g cos 8 dB 

C 

(2.20) 

(2.21) 

B1 = $1 g sin B dB 

Pi 

Then 

(g - E) id0 = 0 6 

and if we take 

it follows that to first order 

&-A gzdB=-X (i+g-i)EdB=-X i2dB<0 6 6 6 

If the flow is subsonic this procedure should converge toward the desired speed distribution 

since the solution will remain smooth, and no unbounded derivatives will appear. If, however, the 

flow is transonic, one must allow for the appearance of shock waves in the trial solutions, even if 

~d is smooth. Then q - qd is not differentiable. This difficulty can be circumvented by a more I 
I 16 

I sophisticated choice of the cost function. Consider the choice 

C 

where X1 and A2 are parameters, and the periodic function S( 0) satisfies the equation 

Then ~ 

fl = J (A1 S6S + x2 ifgay dS &)de 

C 

(2.22) 

(2.23) 

= Sbqd0 

C 

Thus S replaces q - qd in the previous formulas, and if one modifies the boundary condition (2.13) 

to 

the formula for the gradient becomes 

(2.24) 

(2.25) 
17 

i 

1 instead of equation (2.18). Then one modifies f by a step -Ai in the direction of the projected 

gradient as before. 

I 

The final design procedure is thus as follows. Choose an initial profile and corresponding 

mapping function f. Then 

1) Solve the flow equations (2.2 - 2.7) for 4, u, v, q, p. 

2) Solve the ordinary differential equation (2.23) for S. 

i 

I 3) Solve the adjoint equation (2.12) for $ subject to the boundary condition (2.24). 

~ 4) Solve the auxiliary Poisson equation (2.15) for P. 

I 5) Evaluate 

on C, and find its projection i onto the admissible subspace of variations according to 

equations (2.20) and (2.21). 

6) Correct the boundary mapping function fc by I 

I 

~ 

and return to step 1. 
I 18 

I 3. Design for the Euler eauations using conformal maDDing 

This section treats the case of two dimensional compressible flow where the potential flow 

equation is replaced as a mathematical model by the inviscid Euler equations. Let p, p, u, v, E 

and H denote the pressure, density, Cartesian velocity components, total energy and total 

enthalpy. For a perfect gas 

and 

pH = pE + p 

where ris the ratio of specific heats. The Euler equations may then be written as 

where x and y are Cartesian coordinates, t is the time coordinate and 

(3.4) 

As in the previous section, suppose that the domain D exterior to the profile C in the z 

plane is mapped conformally onto the domain exterior to a unit circle in the Q plane (see Figure 

1). Assume also that the outer boundary B of the domain is very far from the profile. Let the 

derivative of the mapping function be 

==he dz iP (3.5) 
19 

9 F= 

Also let r and 0 be polar coordinates in the u plane, where in this case it is more convenient to 

take r as the true radial coordinate denoted by R in the previous section, and 8 is measured in the 

- PU 

puu + sp 

puv -cp 

~ PUH 

clockwise direction. Define the rotation parameters 

s = sin (p- 0) 

and rotated velocity components 

[:]=[:-:I[: 

Then the Euler equations become 

a (rh2w) + a (hF) + 5 a (rhG) = 0 

where 

W= 9 G= 

PV 

pvu + cp 

pvv + sp 

PVH 

(3-7) 

(3.9) 

Then the flow is determined as the steady state solution of equations (3.8) and (3.9), subject to 

the flow tangency condition 

V=O onC (3.10) 

At the far field boundary B conditions can be specified for incoming waves, while outgoing waves 

are determined by the solution. 
20 

In contrast to the case of potential flow, the pressure is not determined solely by the speed, 

and assuming that one wishes to control the surface pressure distribution, a suitable cost function 

is 

(3.11) 

where pd is the desired pressure. A modification to the mapping function Will influence equations 

(3.8) and (3.9) through changes bh and S/3 in both the modulus and argument of 

leading to a variation in the cost function 

, finally 

where 6p is the variation of the pressure. 

Now the mapping variations cause variations in the rotation parameters 

6s = cap, & = - sap 

qhs) = ~6h + hc60, qhc) = c6h - hs6B 

(3.12) 

(3.13) 

Define the Jacobian matrices 

CZSA-CB, D = CA + SB (3.14) 
21 

Then the variation 6w in w satisfies 

dhCbw) a + -aE(rhDbw) a = - dF6h a + hGb@ - & r(G& - hFb0) 

Also 

N=O onC 

(3.15) 

(3.16) 

At the outer boundary there will be no variation in characteristic variables corresponding to 

incoming waves. If we take the outer boundary B at a fixed radius, incoming waves correspond to 

negative eigenvalues of D. Suppose that D is represented as T A T-l, where A is a diagonal 

matrix containing its eigenvalues, and the columns of T are eigenvectors of D. Define 

and 6;- as the components of & corresponding to negative eigenvalues of T. Then 

&-=O onB (3.17) 

Since bw satisfies the constraint (3.15), we can replace equation (3.12) by 

6I = (p -pd)6pdB-l 1 $J T{a (hC6w) + (rhDh)}drdB 

(3.18) -Id $J Ta {%(Fbh + hG6p) + ;r(Gbh-hFbfl}drdB I 

i 

I where the vector $J is a Lagrange multiplier, and the superscript T denotes the transpose. Suppose 
22 

that $ is the steady state solution of the adjoint equation 

2-C T8 $$-rDTg=O inD (3.19) 

At the outer boundary B conditions can be specified for incoming waves, corresponding to positive m 

eigenvalues of DT = T-'" A TT. Define 

+ and '31 as the components of $ corresponding to positive eigenvalues of D. Then we can set 

?+=o on^ (3.20) 

If we integrate equation (3.18) by parts the contribution 

I T rh $ DGwdO = rh? A 66 de 

B 

vanishes because of the complementary boundary conditions (3.17) and (3.20) satisfied by & and 

9 at the outer boundary. If 6h and Gfl decay fast enough in the far field the contribution 

1 1 r (Gbh - hF@) de 

will also be negligible. Thus we find that 

61 = I (p - pd) GpdO + lIT (hD6w + G6h - hF6P) de + J 

I 

I 

C 
23 

where 

J = 11 { (FT$O + GTr$r)6h + (GT$e- FTr$ r ) h@}drdO 

D 

(3.21) 

Also 

Thus using the relations (3.13) 

Now let $ satisfy the boundary condition 

h(c$2 + s$3) = - (P - Pd) on c (3.22) 

Then 

(3.23) 

Finally we can use the fact that the mapping function is fully determined by its boundary 

value to reduce J to a boundary integral. Set 

dz log = f + iP 

f = log l$l = log h 
l and 

Also f and p separately satisfy Laplace’s equation 

~ 

Af=O, A@= 0 

1 and jointly they satisfy the Cauchy Riemann equations 

fg = rPr 7 Po = - rfr 

Let the auxiliary function P satisfy the equation 

AP = h(FTq0 + GTr7CI,) in D 

and the boundary condition 

I P=O onC 

Also let the auxiliary function Q satisfy the equation 

AQ = h(GTtlg - FTrq) in D 

and the boundary condition 

(3.24) 

(3.25) 

(3.26) 

(3.27) 
25 

Then 

J = 51 (AP 6f+ AQ60) drd8 

D 

Thus finally 

SI = 5 g6fd8 

C 

where 

(3.28) 

(3.29) 

As in the previous section, an appropriate modification off is 

where is the projection of g onto the admissible subspace of variations defined by equations 

I (2.20) and (2.21), and X is a sufficiently small positive number. Then 

&=-A g dO<O 1 I5 N2 

C 

I 

If the €low is transonic, shock waves are likely to be formed, and again it may be desirable 

to use a more sophisticated cost function to produce a smooth shape change. In this case we can 

I 
set 

26 

(3.30) 

where XI and A2 are positive parameters, and the periodic function S( 0) satisfies the equation 

Then 

(3.31) 

= J- Sbde 

C 

Thus S replaces p - pd in the previous formulas. If one modifies the boundary condition (3.22) to 

h(cl(12 + sl(13) = - S on C (3.32) 

the formula for the gradient becomes 

g = + 3 - sp 

instead of equation (3.29), and an appropriate modification off is again -Ai. 

(3.33) 
27 

The final design procedure using the Euler equations is thus as follows. Choose an initial 

profile and corresponding mapping function f. Then 

Solve the flow equation (3.8) for w by integrating to a steady state. 

Solve the ordinary differential equation (3.31) for S. 

Solve the adjoint equation (3.19) with the boundary conditions (3.20) and (3.32) for q5 by 

integrating to a steady state. 

Solve the auxiliary Poisson equations (3.24) and (3.26) for P and Q. 

Evaluate 

g = + 3- sp 

on Cy and find its projection g onto the admissible subspace of variations according to 

equations (2.20) and (2.21). 

Correct the boundary mapping function fc by 

6f = -Ai 

where X > 0, and return to step 1. 
28 

4. Wing design using the Euler eauations 

In order to illustrate further the application of control theory to aerodynamic design 

problems, this section treats the case of three-dimensional wing design, again using the inviscid 

Euler equations as the mathematical model for compressible flow. In this case it proves 

convenient to denote the Cartesian coordinates and velocity components by xl, 3, x3 and ul, u2, 

u3, and to use the convention that summation over i = 1 to 3 is implied by a repeated index i. 

The three-dimensional Euler equations may then be written as 

where 

1 and 

Also 

f, = f3 = I pu3 

pu3u 1 

pH = pE + p 

pu3u2 1 

(4.2a) 

(4.2b) 

(4.3) 
29 

Consider a transformation to coordinates X1, X2, X3 where 

axi 1 axi 

J = det(H), Hfj J Hij = 9 9 

The Euler equations can now be written as 

aFi 

at aw+ q=O 

where 

J 

Define the contravariant velocity vector 

u1 

= E-1 Iu2 

Then 

Fi = J 

(4.4) 

(4.5) 

Assume now that the new coordinate system conforms to the wing in such a way that the wing 
surface Bw is represented by X2 = 0. Then the flow is determined as the steady state solution of 

equation (4.5) subject to the flow tangency condition 

U2=0 onBw (4.9) 

At the far field boundary, conditions can be specified for incoming waves as in the 

two-dimensional case, while outgoing waves are determined by the solution. 

Suppose now that it is desired to control the surface pressure by varying the wing shape. It 

is convenient to retain a fixed computational domain. Variations in the shape then result in 

corresponding variations in the mapping derivatives defined by H. Introduce the cost function 

(4.10) 

where Pd is the desired pressure. A variation in the shape will cause a variation Sp in the pressure 

and consequently a variation in the cost function 

Since p depends on w through the equation 01 state (4.3), the variation Sp can be 

determined from the variation h. Define the Jacobian matrices 

ati 

A==, Ci = H..A 1 1J J 

(4.11) 

(4.12) 

Then 

(4.13) 
31 

where 

ax. 

I 

I 

I 

I 

I 

I 

i 

I 

I 

i 

6Fi = CisW + b(J &) f. 

jJ 

and for any differentiable vector $ 

L$6Fidv=jni) T bFids 

1 

bouda r ies 

(4.14) 

(4.15) 

where nl, n2 and n3 are the components of a unit vector normal to the boundary. On the wing 

= 0 and it follows from equation (4.9) that surface Bw , n1 = "3 

Suppose now that $ is the steady state solution of the adjoint equation 

g-Ci&=O Ta inD 

(4.16) 

(4.17) 

At the outer boundary D incoming characteristics for ) correspond to outgoing characteristics for 

6w. Consequently, as in the two-dimensional case, one caa choose boundary conditions for $ such 
32 

that 

rn 

If the coordinate transformation is such that S(JH-') is negligible in the far field, the only 

remaining boundary term is 

Let + satisfy the boundary condition 

(4.18) 

Then, since it follows from equation (4.17) that 

Ci 6w dV = 0 

we find that 
33 

ITY i 

X 

Figure 2 

A convenient way to treat a wing is to introduce sheared parabolic coordinates through the 

transformation 

x = g {x2 - (Y + S(X,Z))2} 

y = X(Y + S(X,Z)) 

z=z 

Here x, y, z are Cartesian coordinates, and X and Y+S correspond to parabolic coordinates 

generated by the mapping 

at a fixed span station Z. The surface Y=O is a shallow bump corresponding to the wing surface, 

with a height S(X,Z) determined by the equation 

X + is = 2 (xs + iy,) 1 

where xS(z) and y,(z) are coordinates of points lying in the wing surface. We now treat S(X,Z) as 
34 

the control. 

In this case 

H= 

while 

X-(Y+S)Sx 

Y+S + xsx 

0 

J = X2 + (Y+S)2 

and 

X 

H-'= $ [ - (Y+S+XSx) 

0 

Also 

6J = 2(Y+S) 6s 

0 

- ( sS+X6SX) 

0 

xsz 1 I 

0 

-JSz 

J 

0 

-(sSSX+(Y+S)6Sx) - (6J Sz+JsSz) 

0 1 
35 

I 

Inserting these formulas in equation (4.19) we find that the volume integral in 61 is 

$$ { ( 6S+X6Sx)fl + ( 6SSx+(Y+S)6Sx)f, + ( 6JSz+JaS,r,)] dV -d 

I I where S and 6s are independent of Y. Therefore, integrating over Y, the variation of the cost 

I function can be reduced to a surface integral of the form 

I 

I 

61 = (P(~,z)& + ~(X,Z)&, + R(X,Z)6SZ)dXdZ 

I Also the shape change will be confined to a bounded region of the X-Z plane, so we can integrate 

by parts to obtain 

a= (P -%-E) 6s dX dZ 

! Thus to reduce I we can choose 

In order to impose a thickness constraint we can define a baseline surface So(X,Z) below 

> which S(X,Z) is not allowed to fall. Now if we take X = X(X,Z) as a non-negative function such 
36 ' that 

S(X,Z) + WGq 2 SO(X,Z) 

Then the constraint is satisfied, while 
37 

5. Conclusion 

The purpose of the last three sections is to demonstrate by representative examples that 

control theory can be used to formulate computationally feasible procedures for aerodynamic 

design. The cost of each iteration is of the same order as two flow solutions, since the adjoint 

equation is of comparable complexity to the flow equation, and the remaining auxjliary equations 

could be solved quite inexpensively. Provided, therefore, that one can afford the cost of a 

moderate number of flow solutions, procedures of this type can be used to derive improved 

designs. The approach is quite general, not limited to particular choices of the coordinate 

transformation or cost function, which might in fact contain measures of other criteria of 

performance such as lift and drag. For the sake of simplicity certain complicating factors, such as 

the need to include a special term in the mapping function to generate a corner at the trailing 

edge, have been suppressed from the present analysis. Also it remains to explore the numerical 

implementation of the design procedures proposed in this paper. 
1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10 

11. 

12. 

13. 

14. 

References 

PRECEDING PAGE BLANK NOT FILMED 

39 

Murman, E.M. and Cole, J.D., "Calculation of Plane Steady Transonic Flows", AIAA 

Journal, Vol. 9, 1971, pp. 114-121. 

Jameson, A., "Iterative Solution of Transonic Flows Over Airfoils and Wings, Including 

Flows at Mach l", Comm. Pure. Appl. Math, Vol. 27, 1974, pp. 283-309. 

Jameson, A. and Caughey, D.A., "A Finite Volume Method for Transonic Potential Flow 

Calculations", Proc. AIAA 3rd Computational Fluid Dynamics Conference, Albuquerque, 

1977, pp. 35-54. 

Bristeau, M.O., Pironneau, O., Glowinski, R., Periaw, J., Perrier, P., and Poirier, G., 

"On the Numerical Solution of Nonlinear Problems in Fluid Dynamics by Least Squares 

and Finite Element Methods (11). Application to Transonic Flow Simulations", Proc. 3rd 

International Conference on Finite Elements in Nonlinear Mechanics, FENOMECH 84, 

Stuttgart, 1984, edited by J. St. Doltsinis, North Holland, 1985, pp. 363-394. 

Jameson, A., Schmidt, W., and Turkel, E., "Numerical Solution of the Euler Equations by 

Finite Volume Methods Using Runge-Kutta Time Stepping Schemes", AIAA Paper 

81-1259, AIAA 14th Fluid Dynamics and Plasma Dynamics Conference, Palo, Alto, 1981. 

Ni, Ron Ho., "A Multiple Grid Scheme for Solving the Euler Equations", AIAA Journal, 

Vol. 20, 1982, pp. 1565-1571. 

Pulliam, T.H. and Steger, J.L., "Recent Improvements in Efficiency, Accuracy and 

Convergence for Implicit Approximate Factorization Algorithms", AIAA Paper 85-0360, 

AIAA 23rd Aerospace Sciences Meeting, Reno, January 1985. 

MacCormack, R.W., "Current Status of Numerical Solutions of the Navier-Stokes 

Equations", AIA A Paper 85-0032, AIAA 23rd Aerospace Sciences Meeting, Reno, January 

1985. 

Jameson, A., Baker, T.J. and Weatherill, N.P., "Calculation of Inviscid Transonic Flow 

Over a Complete Aircraft", AIAA Paper 86-0103, AIAA 24th Aerospace Sciences Meeting, 

Reno, January 1986. 

Jameson, A., "Successes and Challenges in Computational Aerodynamics", AIAA Paper 

87-1184-CP, 8th Computational Fluid Dynamics Conference, Hawaii, 1987. 

Lighthill, M.J., "A New Method of Two Dimensional Aerodynamic Design", ARC, Rand M 

2112, 1945. 

McFadden, G.B., "An Artificial Viscosity Method for the Design of Supercritical Airfoils", 

New York University Report COO-3077-158, 1979. 

Taverna, F., "Advanced Airfoil Design for General Aviation Propellers", AIAA Paper 

83-1791, 1983. 

Garabedian, P. and McFadden, G., "Computational Fluid Dynamics of Airfoils and 

Wings", Proc. of Symposium on Transonic, Shock, and Multidimensional Flows, Madison, 

1981, Meyer, R., ed., Academic Press, New York, 1982, pp. 1-16. 
~ 15. i 

16. 

17. 

18. 

19. 

20. 

21. 

40 

Tranen, J.L., "A Rapid Computer --ided Transonic .A 

74-501, 1974. 

iil Design Method, A 4A Paper 

Volpe, G. and Melnik, R.E., "The Design of Transonic Aerofoils by a Well Posed Inverse 

Method", Int. J. Numerical Methods in Engineering, Vol. 22, 1986, pp. 341-361. 

Henne, P.A., "An Inverse Transonic Wing Design Method", AIAA Paper 80-0330, 1980. 

Garabedian, P.R. and Korn, D.G., "Numerical Design of Transonic Airfoils", Proc. 

SYNSPADE 1970, Hubbard, B., ed., Academic Press, New York, 1971, pp. 253-271. 

Giles, M., Drela, M. and Thompkins, W.T., "Newton Solution of Direct and Inverse 

Transonic Euler Equations", AIAA Paper 85-1530, Proc. AIAA 7th Computational Fluid 

Dynamics Conference, Cincinnati, 1985, pp. 394402. 

Hicks, R.M. and Henne, P.A., "Wing Design by Numerical Optimization", AIAA Paper 

79-0080, 1979. 

Lions, Jacques Louis, "Optimal Control of Systems Governed by Partial Differential 

Equations", translated by S.K. Mitter, Springer Verlag, New York, 1971. 
1. Report No. 2. Government Accession No. 

NASA CR-181749 

ICASE Report No. 88-64 

4. Title and Subtitle 

AERODYNAMIC DESIGN VIA CONTROL THEORY 

3. Recipient's Catalog No. 

5. Report Date 

7. Author(s) 

Ant ony Jameson 

10. Work Unit No. 1 

8. Performing Organization Report No. 

88-64 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

and Engineering 

Mail Stop 132C, NASA Langley Research Center 

Hampton, VA 23665-5225 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 I 

15. Supplementary Notes 

505-90-21-01 

11. Contract or Grant No. 

NAS1-18107 

- 

13. Type of Report and Period Covered 

Langley Technical Monitor: 

Richard W. Barnwell 

17. Key Words (Suggested by Author(s)) 

control theory, aerodynamic design 

Final ReDort 

18. Distribution Statement 

02 - Aerodynamics 

Unclassified - unlimited 

Submitted to J. Sci. Comput. 

19. Security Classif. (of this report) 20. Security Classif (of this page) 21 No. of pages 

Unclassified Unclassified 42 

16. Abstract 

This paper addresses the question of how to modify in aerodynamic design to 

improve the performance. Representative examples are given to demonstrate the 

computational feasibility of using control theory for such a purpose. 

An introduction and historical survey of the subject is included. 

22 Price 

A0 3