The Present Status, Challenges, and Future Developments in Computational Fluid Dynamics Antony Jameson Department of Mechanical and Aerospace Engineering Princeton University Princeton, New Jersey 08544 USA
1. SUMMARY
This paper presents a perspective on corn utational fluid t addresses the dynamics as a tool for aircraft desi requirements for effective industriaf%e, .and tradeoffs between modelling accuracy and computational costs. Issues in algorithmaesi ari discussed in detail, together with a unified,ap roac to the design of shock captunng algonthms. Fin aYly, the pa r discusses the use of techniouec drawn from controf%eow, to detemune ootimal ... aerodynamic shapes. In the future mul$disclplinaj analysis and optimization should be combined to provide an integratednumerical design environment.
P
8"
~~
~~~~~~~
Des ite the advances that have been made, CFD is still not geing exploited as effectively as one would like in the design process. This is p a d due to the long setup and high costs, both human andlcomput?tional of complex flow simulations. The essential requuements for industrial use are: 1. assured accuracy
~~~~~~
~~
~
2. INTRODUCTION Computational methods first began to have a significant impact on aerodynamic analysis and design in the period of 196575. This decade saw the introduction of anel methods which could solve the linear flow mode s for arbitrarily complex geometry in both subsomc and supersonic flow 158, 147, 1791. It also saw the a pearance of the first satisfactory methods for treating &e noAnear equationsof transonic flow [ 123,122,63,64,43,54 and the development of the hodograph method for the dgsign of shock free supercriticalairfoils [U]. Computational Fluid Dynamics (CFD)has now matured to the point at which it is widely accepted as a key tool for aerodynamic desi n Algorithms have been the sub'ect of intensive deviopment for the pa$ two decades. h e principles under1 in the design and mplementauon of robust schemes wLcf~can accurately resolve shock waves and contact discontinuiues in compressible flows are now quite well established. It is also uite well understood how to design high order schemes?or viscous flow, including compact schemes and spectral methods. Adaptive refinement of the mesh interval h) and the order of approximations (p) has been success ully exploited both separately and in combination in the hp method [1261. A continuing obstacle to the treatment of configurations with complex geometry has been the problem of mesh generation. Several eneral techniques have been developed, includin alg&aic transformations and methods based on the sofution of elliptic apd hyperbolic uations. In the last few years methods using unstructure meshes have also be un to gain more general acceptance. The Dassaultled the way in developing a finite element me od or transom potenual flow. They obtained a solution for a complete Falcon 50 as early as 1982 [Z]. Euler methods for unstructured meshes have been the subject of intensive development by several oups since 1985 []IO, 82, 81, 163, 141, and Naviergokes methods on unstructured meshes have also been demonstrated[ll7, 118.111.
P
z
&r~
7
2. acceptable computational and human costs 3. fast turn around.
Improvements are still needed in all three areas. In particular, the fidelity of modelling of high Reynolds number viscous flows continues to be limited by computational costs. Consequently accurate and costeffective simulation of viscous flow at Reynolds numbers associated with full scale Ai ht, such as the prediction of high lift devices, remains a chlenge. Several routes are available toward the reduction of computational costs, includ.in the reduchon of mesh requlrements by the use of hgfer order schemes, improved convergence to a steady state b so histicated acceleration methods, fast inversion me&od; !or im licit schemes, and the exploitation of massively p a r a d computers. Another factor limiting the effective use of CFD is the lack of good interfaces to computer aided design (CAD) systems. The geometry models provided by existing CAD systems oftenfail to meet the riquiremenis of continuity and smoothness needed for flow simulation, with the consequence that they must be modified before they can be used to rovidethe inputfor mesh generation. This bottleneck, w%ch impedes the automation of the mesh generation rocess, needs to be eliplinated, "d the CFD softyare shoufd be fuUy integrated m a numencal desi envuonment. In addiuon to more accurate and coste%ctive flow prediction methods, better optimizations methods are also needed, so that not only can designs be rapid1 evaluated, but duections of improvement can be identiJed. Possession of techniques which,result in a faster design cycle gives a crucial advantage in a competitive envuonment. A critical issue, examined in the next section, is the choice of mathemaucal models. What level of complexity is needed to provide sufficient accuracy for aerodynamic desi n and what is the im act on cdst and turnaround time. %Secuon , . 3 addresses &e design of numerical algorithms for flow simulation. Section 4 resents the results of some numerical calculations whicl require moderate computer resources and could be completed with the fast turnaround required by industrial users. Section 5 discusses automatic desien omedures which can be used to produce opumum a&dynamic designs. Finally, Secnon 7. offers an outlook for the future.
Paper presented ut the AGARD FDP Symposium on "Progress and Challenges in CFD Methods and Algorithms" held in Seville, Spuin, from 25 October 1995, and published in CP578.
12
3. THE COMPLEXITY OF FLUID FLOW AND MATHEMATICAL MODELLING 3.1 The Hierarchy of Mathematical Models Many critical henomena of fluid flow, such as shock waves and turI? ulence,, are essentlally nonlinear. They also exhibit extreme disparities of scales. m l e the actual thickness of a shock wave is of the order of a me? free path of the gas articles, on a macrosco ic scale its thickness . i. s essenti$Iv, zero. In turbulent ow enerev is transferred from large scale motions to pro essive? sm4ler eddies until the scale becomes so sma I that the motlon motion is dissipated by viscosity. The ratlo of the length perbistlng scale of the global flow to that of the smallest persisting eddies is of the order Rei, where Re is the Reynolds numher, typically in the range of 30 million for an aircraft. In order to resolve such scales in all three space directions a computational grid with the order of Rei cells would be required. This is beyond the ran e of any current or foreseeable computer. Consequeniy mathematical models with varyin degrees of simplification have to be introduced in ode, to make corn utational simulatlon of flow feasible, and to produce viagle and costeffective meth~~
~~~~
~~~
~~~~~
~
l
f
ods.
Figure I (su plied by Pradee Ra') indicates a hierarchy of mode[ at different leve% o(simp1ification which have Droved useful in practice. Efficient flight is generally achieved by theuse of smooth and
potential flow or Euler solutions for an airfoil can be accurately calculated on a mesh with 160 cells around the section, and 32 cells normal to the section. Using multigrid techniques 10 to 25 cycles are enou h to Gbtain a conver ed result. Consequently airfoil c&uIations can be perkrmed in seconds on a Cray W,and can also be performed on 486classdpersonal computers. Correspondingl accurate three mensiond mviscid calculanons can ge performed for a wing on a mesh, say with 192x32~48=294,912 cells, in about 5 minutes on a single processor Cra YMP, or less than a minute with eight processors, or in or 2 hours on a workstation such as a Hewlett Packard 735 or an IElM 560 model. VISCOUS simulations at high Reynolds numbers require Careful twodimensional studies
1
ary layer, in addition to 32 intervals between the boundary layer and the far field, leading to a total of 64 intervals. In order to prevent degradations in accuracy and converence due to excessively large aspect ratios (in excess of f.000) in the surface mesh cells, +e chordwise resolution must also be increased to 512 intervals. Reasonably accurate solutions can be obtained in a 512x64 mesh in 100multigrid cycles. Translated to three dimensions, this would imply the need for meshes with 510 million cells (for example, 512x64x256= 8,388,608 cells as shown in Figure 2). When simulations are performed on less fine meshes with, sa ,500,000 to 1million cells, it is very hard to avoid mesh Jependency in the solutionsas well as sensitivity to the turbulence model.
p.
t
Figure 2 Mesh Requirements for a Viscous Simulation
Figure 1: Hierarchy of Fluid Flow Models
3.2 Computational Costs Computational costs vary drasticall with the choice of mathematical model. Panel methois can t x effectively used to solve the linear potential flow equation with higherend personal com uters (with an Intel 80486 micro rocessor, for exampL). studies of the dependency of %e result on mesh refinement, performed by this author and others, have demonstrated that inviscid transonic
A typical algorithm requires of the order of 5,000 floating oint o ratlonsper mesh point in one multigrid iteration. b t h l!$milllion mesh points, the operationcount is of the order of 0 . 5 10" ~ per cycle. Given a computer capable of sustaining IOii operationsper second (100 gi aflo s), 200 cycles could then be performed in 100 secon%. &Uulations of unsteady viscous flows (flutter, buffet) would be likely to re uire l,OO010,000 tune steps. A further progression to arge eddy simulation of complex c o d uratlons would require even eater resources. The f$ lowing estimate is due to W.# Jou [90].Suppose that a conservative estimate of the size of edkes in a boundary layer that ought to be resolved is 1/ 5 of the boundary layer thickness. Assuming that IO points are needed to resolve
4
13
a single eddy, the mesh interval should then be 1/50 of the boundary layer thickness. Moreover, since the eddies are threedimensional, the same mesh interval should be used in all three directions. Now, if the boundary layer thickness is of the order of 0.01 of the chord length, 5,000 intervals will be needed in the chordwise direction, and fora win with an aspect ratjo of IO,50,000 intervals will be needA in the spanwise direction. Thus, of the order of 50 x 5 , 000 x 50, 000 or 12.5 billion mesh points would be needed in the boundary layer. If the time dependent behavior of the eddies is to be fully resolved using time steps on the order of the time for a wave to pass through a mesh interval, and one allows for a tolal lime e ual to the time required for waves to travel three times ?he length of the chord. of the order of 15,ooO time steps would be needed. Performance beyond the teratlop ( IOi2 operations per second) will he needed to attempt calculations of this nature, which also have an informauon content far beyond what is needed for enginering anal sis and design. The desi er does not need to know x e details of the eddies in boundary layer. The primary purpose of such calculations is to im rove thc prediction of avera ed uantities such as skin kction, and the prediction of gfobagbehaviorsuch as the onset of separation. The man current use of NavierStokes and large eddy simulations is to gain an improved insight into ~e ph sics of turbulent flow, which may in tum lead to the deve opment of more comprehensive and reliable turbulence models.
tE
r
The selection of sufficiently accurate mathematical models and a jud ent of their costeffectiveness ultimately rests with i n g q Aircraft and spacecraft desi ns nor mally pass throug the three phases of conce tuafdesigi prelimnary design, and detailed desi n torrespond!ngl, the appropriateCFD models witviry in complexity. the Conceptual and prelinary desi n phases, the emphasis will be on relatively smple which can give results with very rapid turnaround and low computer costs. in order to evaluate alternative confieurationc and
&
modis
Ild be placed onnumerical p r e c tions bas forced the extensive use of wind tunnel tesung at an early stage of the desi n This practice was very exoensive. The limited n u d e ; of models that could be fahcated also limited the range of design variations that could be evaluated. It can be anticipated tha! in the future, the role of wind tunnel tesung in the design process will be more one of verification. Experimena &search to improve our understanding of the physics of complex flows will continue, however, to play a vital role. ~~~~~
~
~~~
CFDALGORlTHMS 4.1 Difficulties of Flow Simulation 4.
The corn utational simulation of fluid flow presents a number of' severe challenges for al orithm design. At the level of inviscid modeling, the in%erent nonlinearity of the fluid flow equations leads to the formation of singularities such as shock waves and contact discontinuihes. Moreover, the geometric configurations of interest are extremely com lex, and generally contain sh edges which lead to %e shedding of vortex sheets. ?xtreme ients near stagnation points or win tips may also to numerical errors that can have iobal influence. Numericallv generated entrow mav be convected from the leadingedge, for cxampl'e: cauiing the formauon of a numericall induced boundary layer which can lead to separahon. &e need to treat extenor domans of infinite extent IS also a source of difficulty. Boundary condihons imposed at artificial outer boundahes may cause reflected waves which significantly interfere with the flow. When viscous effects are also included in the simulation, the extreme difference of the scales in the viscous boundary layer and the outer flow. which is essential1 inviscid, IS another source,ofdifficulty. forcing the use o meshes with extreme variations in the mesh intervals. For these reasons, CFD has been a driving force for the development of numerical algorithms.
3.3 'Aubnlence Modelling It is doubtful whether a universally valid turbulence model, capable of describing all complex flows, could he devised [52]. AI ebraic models [30,9] haveprovedfairly satisfactory for %e calculation of attached and slightly separated wing flows. ,These models rely on the boundary layer concept, usual1 incorporating separate formulas for the m e r and outer a ers, and they require an estimate of a length scale whici depends on the thickness of the boundary layer. The estimation of this quantity by a search for a maximum of the vorticity times a distance to the wall, as @ the BaldwinLomax model, can lead to ambiguities in internal flows, and also m complex vorhcal flows over slender bodies and highly swept or delta wings 140, 1151. The JohnsonKin model [88], which allows for nonequilibrium effects &rough the inaoduction of an ordinary differential equation for the maxmum shear stress. has improved the rediction of flows with shock induced separation [148, $11, Closure models depending on the solution of transport eouations are widelv acceoted for industrial anolications. f i e s e rnodelseliminatethk need toesumate a llngescale by detecting the ed e of the bounday layer Eddy viscosity models typical y use two equauons for the turbulent 4.2 Structurpd and Unstructured Meshes kinetic energy k and the dissipation rate E , or a pair of The al oritbm designer faces a number of critical deciuivalent quantities [89, 178. 160, 1, 121, 351. Models sions. %%efirst choice that must be made IS the nature !his type enerall tend to resent difficulties in the of the mesh used to divide the flow field into discrete region very cfose to &e wall. &ey also tend to be badly subdomains. The discretization procedure must allow for conditionedfornumerical solution. The ICI model [154] the treatment of complex configurations. The principal is designed to alleviate this problem by taking advantage alternatives are Cartesian meshes, bodyfitted curvilinear of the linear behaviour of the leneth scale 6 near the wall. meshes, and unstructured tetrahedral meshes. Each of In an alternative a proach to theldesi of models which these a proaches has advantages which have led to their are more amenahre to numerical sogtion, new models use. d e Cartesian mesh minimizes the corn lexity of requiring the solution of one transport equation have rethe algorithm at interior points and facilitates %e use of centl been introduced [IO, 1591. %e performance of high order discretization procedures, at the expense of the alygebraic models remains competlhve for wing flows, greater corn lemy, and possibly a loss of accuracy, in the but the one and twoequation models show promise for treatment oPboundary,conditions at curved surfaces. his broader classes of flows. In order to achieve greater unidifficultymay be alleviated b using mesh refinement ro versality, research is also bein ursued on more complex oedures near the surface. d h theu aid, schemes w f i c i Reynolds stress transport m&i, which require the soluuse Cartesian meshes have recentl been develo ed to tion of a larger number of transport equations. treat very complex configurations 149,22,9!]. Another direction of research is the attempt to devise have been widely used and are parmore rational models via renormalization group (RNG) to the treatment of viscous flow betheory [182,155 Both algebraic and twoequation k allow the mesh to be compressed near models devised by this approach have shown promising W~ this approach, the problem of results [I 161. mesh generation itself bas proved to be a major pacing
7
e
r
P
8
rho,
14
equations. In the finite volume method [llZd, the. discretization is accomplished by dividing the omam of the flow into a large number of small subdomans. and ~170,171,157,158],andme~odsbasedonthesolutionof applying the conservation laws in the integral form yperbolic equahons marehmg out from,the body [161 In order to treat very complex codguratlons it general y ent to use a multiblock 1177, 150 roce proves dure, wi separately generated meshes m eact! {loci which may then be atcbed at block faces, or allowed Here f is the flux appearing in equation (1) and dS is to overlap, as m the ghimera scheme [19,20]: While a the directed surface element of the boundary aR of the number of interactive software systems for d eneradomain R. The use of the integral form has the advantage tion have been developed, such as EAOLE, %UbGEN. that no assum tion of the differentiability of the solutions and ICEM, the generahon of a sapsfacto grid for a very is implied, wik the result that it remains a valid statement complex configuratlon may requue mon s of effort. for a subdomain containinga shock wave. In general the The alternative is to use an unstructured mesh in which the subdomains could be arbitrary, but it is convenient to use domainis subdwidedinto tetrahedra. Thisin turnrequlres either hexahedral cells m a body confornun curvilinear the development of solution algorithms ca able of yield mesh or tetrahedrons 111 an unstructured mes%. ing the required accupcy on unstructure8meshes. This approach has been gaming acceptance, as it is becoming Alternative discretization schemes may be obtained by storin flow variables at either the cell centers or the vera parent that it can lead to a speedup and reduction in tices. h e s e variations are illustrated in Figure 3 for the $e cost of mesh generation that more than offsets the increased complexlty and cost of the flow simulations. 30 twodimensional case. With a cellcentered scheme the discrete conservation law Iakes the form competing procedures for eneratingtnangulahons whch l 1 are Delaunay triangulahon have both proved successi [41, 111, based on concepts introduced at the beginning of the cenhuy by Voronoi [175], and the moving front method [ 1111. item. The most commonly used procedures are algebraic transformations [7,44.46,156], methods based on the solution of elliphc equations, pioneered by Thompson
2.
""P
x
4.3 Finite Difterence, Finite Volume, and Finite Element sehemes Associated with choice of mesh m e is the formulation of thediscretizationprocedurefor~~equationsoffluidEow, which can be expressed as differential conservation laws. In the Cartesian tensor notation. let z,be the coordinates. p, 0, T,and E the pressure, density, temperature, and iod energy, and ut d e velocity corn onents Usin the convention that summation over i=1.13 is imlied%y a repeated subscript j , each ConseiVation equati6n has ihe form ~~~~
~~~~
aw + ='af. o. at
3a: Cell Centered Scheme.
(1)
82,.
For the mass equation
3b: Vertex Scheme.
w=p, fj=puj.
Figure 3 Structured and Unstructured Discretizations.
For the i momentum equation d wv+ dt
where U,, is the viscous stress tensor. For the energy equation w=pE, f 3 = ( p E+p) U,
aT
 U3kUk  l C z v 3
where nis the coefficient of heat conduction. The pressure is related to the density and energy by the equation of state 1
P=(r 1) p (E  p U , )
(2)
in which y is the ratio of specific heats. In the NavierStokes equations the viscous stresses are assumed to be linearly proportionalto h e rate of strain, or
where p and X are the coefficients of viscosity and bulk viscosity, and usually X=2p/3. The finite difference method, which requires the use of a *esian or a spuctured curvhear mesh, directly appromates the hfferenhal operators appearing in these
f.S=O,
(4)
faces
where V is the cell volume, and f is now a numerical estimate of the flux vector through each face. f may be evaluated from values of the flow variables in the cells separated by each face, using upwind biasing to allow for the+ections,of wavepro agahon. With.hexahedral cells, equatlon (4) is very S&I to a finite Merence scheme in curvilinear coordinates. Under a transformation to curvilinear coordinates (j. equation (1) becomes (5)
[e].
where J is the Jacobian determinantof the transformat&on The transformed flux J$ f, corresponds matrix to the dot product of the flux f with a vector face area while J represents the transformation of the cell volume. The finite volume form (4) has the advantages that it is valid for both structuredand unstructured meshes, and that it assuresthat a umfonn flow exactly satisfies the equations, because CfWsS=Ofor a closed control volume. Finite dif€erence schemes do not necessarily satisfy
Jg,
this constraint because of the discretization errors in evaluating and the inversion of the transformation matrix. zj A cellvertex finite volume scheme can be derived by taking the union of the cells surroundinga given vertex as the control volume for that vertex [55, 71, 1391. In equation (4), V is now the sum of the volumes of the surrounding cells, while the flux balance is evaluated over the outer faces of the polyhedral control volume. In the absence of upwind biasing the flux vector is evaluated by averagin over the corners of each face. This has the advantage of remaining accurate on an irregular or unstructured mesh. An alternative route to the discrete equations is provided by the finite element method. Whereas the finite difference and finite volume methods a proximate the differentia1 and integral operators, the nite element method proceeds by inserting an approximate solution into the exact equations. On multiplying by a test function q5 and integratingby parts over space, one obtains the weak form
maximum V k  v j 5 0, and at a minimum V k  v j 2 0. Thus the condition cjk
4.4 Nonoscillatory Shock Capturing Schemes 4.4.1 Local Extremum Diminishing (LED) Schemes
The discretization rocedures which have been described in the last section read to nondissipative approximations to the Euler equations. Dissipative terms may be needed for two reasons. The first is the possibilit of undamped oscillatory modes. The second reason is txe need for the clean capture of shock waves and contact discontinuities without undesirable oscillations. An extreme overshoot could result in a negative value of an inherently positive quantity such as the pressure or density. The next sections summarize a unified ap roach to the constructlon of nonoscillatory schemes via tie introduction of controlled diffusive and antidiffusiveterms. This is the line adhered to in the author's own work. The development of nonoscillato schemes has been a rime focus of algorithm researchxr com ressible flow. Eonsider a general semidiscrete scheme o the form
P
(7) A maximum cannot increase and a minimum cannot decrease if the coefficients C j k are nonnegative, since at a
(8)
is sufficient to ensure stability in the maximum norm. Moreover, if the scheme has a compact stencil, so that c j k ' o when j and k are not nearest neighbors, a local maximum cannot increase and local minimum cannot decrease. This local extremum diminishing (LED) roperty prevents the birth and growth of oscillations. ' h e onedimensional conservation law
au a
R
which is also valid in the presence of discontinuitiesin the flow. In the Galerkin method the a proximate solution is expandedin terms of the same famify of functions as those from which the test functions are drawn. By choosing test functions with local su port, separate equations are obtained for each node. €o!r example, if a tetrahedral mesh is used, and q5 is piecewise linear, with a nonzero value only at a sin le node, the equations at each node have a stencil whick contains only the nearest neighbors. In this case the finite element approximation corres onds closely to a finite volume scheme. If a iecewise Tinear approximation to the flux f is used in tie evaluation of the integrals on the right hand side of equation (6), these integrals reduce to formulas which are identical to the flux balance of the finite volume scheme. Thus the finite difference and finite volume methods lead to essential1 similar schemes on structured meshes, while the finite voyume method is essentially equivalent to a finite element method with linear elements when a tetrahedral mesh is used. Provided that the flow equations are expressed in the conservation law form (l), all three methods lead to an exact cancellation of the fluxes through interior cell boundaries, so that the conservative ro erty of the equations is preserved. The important r o l orthis propert in ensuring correct shock jump conditions was pointedbut by Lax and Wendroff [97].
2 0, k Pj
=o
 + f(U) at ax
provides a useful model for analysis. In this case waves are propagated with a speed a(u> and the solution is constant along the characteristics %=a(u). Thus the LED property is satisfied. In fact the total variation
=%,
of a solution of this equation does not increase, provided that any discontinuit appearing in the solution satisfies an entropy condition [ 61. Harten pro osed that difference schemes ought to be designed so tiat the discrete total variation cannot increase [56]. If the end values are fixed, the total variation can be expressed as
B
T V ( ~=2)
(E maxima  E minima) .
Thus a LED scheme is also total variation diminishing (TVD). Positivit conditions of the ty e expressed in equations (7) a n i (8) lead to diagonaiy dominant schemes, and are the key to the elimination of improper oscillations. The positivity conditions may be realized by the introduction of diffusive terms or by the use of upwind biasing in the discrete scheme. Unfortunately, they may also lead to severe restrictions on accuracy unless the coefficients have a complex nonlinear dependence on the solution. 4.4.2 Artijicial Difision and Upwinding
Following the pioneering work of Godunov [ 5 13, a variety of dissipative and upwind schemes designed to have good shock capturing properties have been developed during the past two decades [162, 23, 98, 100, 146, 130, 56, 129, 166, 5, 68, 183, 62, 180, 13, 12, 111. If the onedimensional scalar conservation law
av a
=o
 + f(v) at ax
(9)
is represented by a three point scheme dvjC.+t
dt
+
(Vj+,
 Vj) + cT
3
(U+,
 Vj) ,
34
the scheme is LED if
c3+.+ t >  0,
2 0.
A conservative semidiscrete approximation to the onedimensional Conservation law can be derived b subdividing the line into cells. Then the evolution of tie value v j in the jth cell is given by
dvj Ax dt + hj+f  hjf=O,
16
where hj+l is an estimate of the flux between cells j and + 1. Th2esimplest estimate is the arithmetic average (fj+l + fj) /2, but this leads to a scheme that does not satisfy the positivity conditions. To correct this, one may add a dissipative term and set
j
In order to estimate the required value of the coefficient aj+4, let aj+4 be a numerical estimate of the wave speed Af
4.4.3 Hi h Resolution Switched Schemes: JamesonScfmidtTurkel (JST) Scheme
Higher order nonoscillatory schemes can be derived by introducing antidiffusive terms in a controlled manner. An early attem t to roduce a hi h resolution scheme by this approac\ is tRe JamesonlchmidtTurkel (JST) scheme [ 8 5 ] . Suppose that antidiffusive terms are introduced by subtracting neighboring differences to produce a third order diffusive flux d 3. + l1= c r3. +2l { A u 3. + I7
 21 (Auj+* + A u j  ~ ) } ,
(15)
&.
which is an approximation to ; a h 3 The positivity condition (8) is violated b this scheme. It proves that it generates substantial oscjlrations in $e vicinity of shock waves, which can be eliminated by switchinglocally to the first order scheme. The JST scheme therefore introduces blended diffusion of the form
Then
where
Avj+f=vj+l  v j 1
and the LED condition (10) is satisfied if
If one takes
The idea is to use variable coefficients and E (4) .+ , 3+4 J T which produce a low level of diffusion in regions where the solution is smooth, but prevent oscillations near discontinuities. If is constructed so that it is of order
1 a.+1=3 I 2 laj+jll
3+4
one obtains the first order upwind scheme
Ax2 where the solution is smooth, while E$; is of order unity, both terms in dj+ 4 will be of order Ax3.
This is the least diffusive first order scheme which satisfies the LED condition. In this sense upwinding is a natural approach to the construction of nonoscillatory schemes. It may be noted that the successful treatment of transonic otential flow also involved the use of U wind biasing. h i s was first introduced by Murman and ??ole to treat the transonic small disturbance equation [ 1231. Another important re uirement of discrete schemes is that they should excluie nonphysical solutions which do not satisfy appropriate entropy conditions [95], which require the convergence of characteristics towards admissible discontinuities. This places more stringent bounds on the minimum level of numerical viscosity [113, 169, 128, 1311. In the case that the numerical flux function is strictly convex, Ais0 has recently proved [2] that it is sufficient that
for e > 0. Thus the numerical viscosity should be rounded out and not allowed to reach zero at a point where the approaches zero. This justifies, for wave speed a ( u ) example, Harten's entropy fix [56]. Higher order schemes can be constructed b introducing higher order diffusive terms. Unfortunate{; these have larger stencils and coefficients of varying si n which are !D scheme, not compatible with the conditions (8) for a E and it is known that schemes which satisfy these conditions are at best first order accurate in the nei hborhood of an extremum. It proves useful in the fol owin development to introduce the conce t of essentially focal extremum diminishing (ELED) scgemes. These are defined to be schemes which satisfy the condition that in the limit as the mesh width Ax + 0, local maxima are nonincreasing, and local minima are nondecreasing.
=%
The JST scheme has proved very effective in practice in numerous calculations of com lex steady flows, and conditions under which it could ge a total variation diminishing (TVD) scheme have been examined by Swanson and Turkel [165]. An alternative statement of sufficient conditions on the coefficients E !'Iand E !4) for the JST 3+t 3+i scheme to be LED is as follows:
Theorem 1 (Positivity of the JST scheme) Suppose that whenever either vj+l or vj is an extremum the coefficients of the JST scheme satisfy 3+t
> 1
"'+f
21
3
I
1
j+io. 
e (4)
Then the JST scheme is local extremum diminishing (LED). Proof: We need only consider the rate of change o v at extrema1 points. Suppose that v j is an extremum. d e n
and the semidiscrete scheme ( 1 1 ) reduces to
B
and each coefficient has the required sign. U
17
In order to construct 63!*) and E 3!4) 4 with the desired prop# erties define
=[
R(u,v)
Ii'%%l
9
if U P O or v #O 0 ifu=v=O,
(18)
Set
and
where q is a positive integer. Then R(u,v ) =1 if U and v have opposite signs. Otherwise R(u,v) < 1 . Now set
Qj=R
Qj+
;=max ( Q j
,Q j + l ) .
Then,
and Thus the scheme satisfies the LED condition if aj+f2 4.4.4 Symmetric Limited Positive (SLIP) Scheme
An alternative route to high resolution without oscillation is to introduce flux limiters to uarantee the satisfaction of the positivity condition (8). e use of limiters dates back to the work of Boris and Book [23]. A particularly sim le wa to introduce limiters, pro osed by the author in 4984 [l8],is to use flux limited issipatton. In this scheme the third order diffusion defined by equation (15) is modified by the insertion of limiters which produce an e uivalent three point scheme with positive coefficients. ?%e original scheme [68]can be improved in the followin manner so that less restrtctive flux limiters are required! Let L ( u , v ) be a limited average of U and v with the following properties:
5%
B
P1. L(U, V I = L ( v ,21)
laj+jI for all j , and $ ( T ) >_ 0, which is assured by property (P4)on L. At the same $me it follows from property (P3)that the first order diffusive flux is canceled when AV is smoothly varying and of constant sign. Schemes constructed by this formulation will be referred to as symmetric limited ositive (SLIP) schemes. This result may be summarize as
B
Theorem 2 (Positivity of the SLIP scheme) Su pose that the discrete conservation law ( I I ) contains a lmited d'ffusiveflux as de ned by equation (20). Then the positivi condition ( I ), together with the prvperties (PIlVYfor limited averages, are sufJicient to ensure satisfaction of the LED principle that a local m i m u m cannot increase and a local minimum cannot decrease. 0
di
A variety of limiters may be defined which meet the requirements of properties (PlP4). Define
P2. L ( a u , av)= a L ( u , v )
P3. L(U, U ) =U P4. L ( U , v ) =O if U and v have opposite signs: otherwise L ( u ,v ) has the same sign as U and U. Properties (PlP3) are natural properties of an average. Pro erty (P4)is needed for the construction of a LED or TVb scheme. It is convenient to introduce the notation $(TI
= L ( l , r )= L ( r ,1 1 ,
1
S ( u , v ) = {sign(u) + sign(v) } 2 which vanishes is U and v have opposite signs. Then two limiters which are appropriate are the following wellknown schemes:
2. Van Leer:
where according to (P4)$ ( T I 2 0. It follows from (P2) on setting a=; or that In order to produce a family of limiters which contains these as special cases it is convenient to set Also it follows on setting v = l and U=T that L(U,V)
Thus, if there exists T < 0 for which $ ( T I > 0, then $ (:) < 0. The only way to ensure that $ ( T I 2 0 is to re uire $ ( T ) =O for all T < 0, corresponding to property (PI). Now one defines the diffusive flux for a scalar conservation law as
1 2
=  D ( u , v ) ( U + w) ,
where D ( U , v ) is a factor which should deflate the arithmetic average, and become zero if U and v have opposite signs. Take
where R(u,v ) is the same function that was introduced in the JST scheme, and q is a positive integer. Then D ( u , v) =O if U and v have opposite signs. Also if q=1, L(u, v) reduces to minmod, while if q=2, L ( u ,v) is
18
equivalent to Van Leer’s li+ter. By generate a se uence of 1imted.averages a limit define by the arithmetic mean when U and v have opposite signs. When the terms are re rouped, it can be seen that with this limiter the SLIP sc eme is exact1 equivalent to the JST scheme, with the switch is definecras
a
Q
one recovers a standard if a j + l < 0. If cr high resolution upwind scheme in semidiscrete form. Consider the case that aj+i > 0 and a j  1 > 0. If one sets + “vj+t  Avj 9 T =T =AV+ 4 ’ A V ~ 4 ’ the scheme reduces to
j+l
€ (4)
=
aj+j (1  Qj+t)
.
d
This formulation thus unifies the JST and SLIP schemes. 4.4.5 Essentially Local Extremum Diminishing (ELED) Scheme with SoftLimiter
The limiters defined by the formula (22) have the disadvantage that they are active at a smooth extrema, reducing the local accuracy of the scheme to first order. In order to prevent.this, the SLIP scheme can be relaxed to give an essentlally local extremum diminishin (ELED) scheme which is second order accurate at smoot extrema by the introduction of a threshold in the limited average. Therefore redefine D (U, U > as
a
D (U,VI =1 
Iq,
(23)
max( IuI + 1v1 ,€Axr>
4,
where T= q 2 2. This reduces to the previous definition if lul + 1vl > €AxT. In any region where the solution is smooth,Avj+s Avj 4 is of order Ax2. In fact if there is a smooth extremum in , series expansion the neighborhood of vj or v ~ +a ~Taylor indicates that Avj+ 3, Avj+ 1 and Avj ;are each individually of order Ax2, since $=O at the extremum. It may be verified that second order accuracy is preserved at a smooth extremum if q 2 2. On the other hand the limiter acts in the usual way if IAvj+ j or lAvj j > €AxT, and it may also be verified that in the limit Ax + 0 local maxima are non increasin and local minima are non decreasing [79]. Thus the sci?eme is essentially local extremum diminishing (ELED). The effect of the “soft limiter” is not only to improve the accuracy: the introduction of a threshold below which extrema of small amplitude are accepted also usually results in a faster rate of convergence to a steady state, and decreases the likelyhood of limit cycles in which the limiter interacts unfavorably with the corrections produced by the updatin scheme. In a scheme recently proposed by Venkatakrisnan a threshold is introduced precisely for this purpose [ 1741.
I
To assure the correct sign to satisf fhe LED criterion the flux limiter must now satisfy the a ditlonal constraint that $ ( T I 5 2. The USLIP formulation is essentially e uivalent to standard upwind schemes [ 130,1661. Both t e SLIP and USLIP constructions can be implemented on unstructured meshes [75,79]. The antidiffusive terms are then calculated by taking the scalar product of the vectors definin an edge with the gradient in the adjacent upstream an downstream cells.
R
9
4.4.7 S stems o Conservation Laws: Flux Splitting and duxD d r e n c e Splitting
Steger and Warming [ 1621first showed how to generalize the concept of upwinding to the system of conservation laws aw a +f(w)=O a t ax by the concept of flux splitting. Suppose that the flux is and have positive and split as f =f + f  where negative ei envalues. Then the first order upwind scheme is producef by taking the numerical flux to be
E
+
This can be expressed in viscosity form as
I
4.4.6 Upstream Limited Positive (USLIP) Schemes
=
 dj+f,
2 where the diffusive flux is
Roe derived the alternative formulation of flux difference s littin [ 1461 by distributing the corrections due to the &x dilference in each interval upwind and downwind to obtain dwj dt
Ax
By adding the antidiffusive correction purely from the upstream side one may derive a family of upstream limited ositive (USLIP) schemes. Corespondin to the original ELIP scheme defined by equatlon (20), akSLIP scheme is obtained by setting
1 (fj+l + f j )
+(fj+l
 f j >  + < f j f j  I >
+ =o,
where now the flux difference fj+l  fj is split. The correspondingdiffusive flux is
dj+f=aj+t {Avj+j  L @vj+f’Avjt)}
if aj+t > 0, or
Following Roe’s derivation, let A j + f be a mean value Jacobian matrix exactly satisfying the condition fj+1
 fj==Aj+i(wj+l  wj) .
(26)
19
Aj+i may be calculated by substituting the weighted averages
into the standard formulas for the Jacobian matrix A=%. A splitting according to characteristic fields is now obtained by decomposing Aj+4 as
Thus these schemes are closely related to schemes which introduce separate s littings of the convective and ressure terms, such ase!lt waveparticle scheme [141,8[ the advection upwind splitting method (AUSM) [ 106, 1761, and the convective upwind and split pressure (CUSP) schemes [76]. In order to examine the shock capturing properties of these various schemes, consider the general case of a first order diffusive flux of the form 1
A ~ +=TAT l ,
(28)
where the columns of T are the eigenvectors of Aj+4, and A is a diagonal matrix of the eigenvalues. Now the corresponding diffusive flux is
where
I A ~ + ~I
=T 1
~ T  11
where the matrix Bj+4determines the properties of the scheme and the scaling factor aj+t is included for convenience. All the previous schemes can be obtained by representing Bj+jas a polynomial in the matrix Aj+t defined by equation (26). Schemes of this class were considered by Van Leer [99]. According to the CayleyHamilton theorem, a matrix satisfies its own characteristic equation. Therefore the third and hi her powers of A can be eliminated,and there is no loss o4enerality in limiting Bj+ to a polynomial of degree 2,
and 1A1 is the diagonal matrix containing the absolute values of the eigenvalues. 4.4.8 Alternative Splittings
Characteristic splitting has the advantages that it introduces the minimum amount of diffusion. to exclude the rowth of local extrema of the charactenstic variables,and Bat with the Roe lineafizatjon it allows a discrete shock structure with a single interior point. To reduce the computational complexity one may replace IAl by a1 where if Q is at least equal to the spectral radius max IX(A) I, then the positivity conditions will still be satisfied. Then the first order scheme simply has the scalar diffusive flux 1
(29)
dj+f=TQj+tAwj+t.
The JST scheme with scalar diffusive flux ca waves with about 3 interior oints, and it has used for transonic flow ca culations because it is bo robust and computationally inexpensive. An intermediate class of schemes can be formulated by defining the first order diffusive flux as a combination of differences of the state and flux vectors
(34)
dj+f=Taj+tBj+f(wj+l  wj) >
Bj+t=aol+ alAj+j + ~ 2 A : + f .
(35)
The characteristic upwind scheme for which Bj+l= substituting Aj+4=TAT', A2 =TA2Tl. Then QO, 3+f ~ 1 and , a2 are determined from the three equations QO
+ Q ~ X ; = I X ~ ~ k=1,2,3. ,
+
The same re resentation remains valid for three dimensional flow gecause Aj+i still has only three distinct eigenvalues U , U + c, U  c. 4.4.9 Analysis of Stationary Discrete Shocks
P
1 .
dj+t=y"j+lc
(wj+l
wj)
1 + z ~ j + 4( f j + l 
fj)
9
(30)
!
I
! WL
i
WL
I!
'vy I
! I
I I where the factor c is included in the first term to make ~ j * + and ~ bj+* dimensionless. Schemes of this class I I are fully upwind in supersonic flow if one takes Q ~ * + ~ = O I I j+l I j+2 I and Pj+t=sign ( M ) when the absolute value of the Mach number M exceeds 1. The flux vector f can be decomFigure 4: Shock structure for single interior point. posed as f = U w + f,, (31) The ideal model of a discrete shock is illustrated in fi where ure (4). Suppose that W L and W R are left and rigat states which satisfy the jump conditions for a stationary f,= (32) shock, and that the corresponding fluxes are fL=f ( W L ) and f ~ (wR). Since the shock is stationary fL= R . . The idea discrete shock has constant states W L to the eft Then and W R to the ri ht, and a single point with an intermediate value W A . e intermediate value is needed to allow f j + i  fj=fi ( w j + l  wj) + tij ( u j + l  u j ) + fpj+,  f p j 7 the discrete solution to correspond to a true solution in (33) which the shock wave does not coincide with an interface where I and U, are the arithmetic averages between two mesh cells. I
( )
f
?a
f
Schemescorrespondingto one, two or three terms in equation (35) are examined in [80]. The analysis of these three
110
cases shows that a discrete shock structure with a single interior point is su ported by artificial diffusion that satisfies the two cond!tions that 1. it produces an upwind flux if !he flow is determined to be supersonic through the interface
the corresponding error in the tem erature may lead to a wrong prediction of associated e ects such as chemical reactions. The source of the error in the stagnation enthalpy is the discrepancy between the convective terms
8
2. it satisfies a generalized ei envalue problem for the exit from the shock of thejorm
(AAR OARBAR) (WR  WA)=O,
(36)
where AAR is the linearized Jacobian matrix and B R is the matrix defining the diffusion for the interface AIR. This follows from the equilibrium condition h R A = h R R for the cell j + 1 in fi ure 4. These two conditions are satisfied by both the ciaracteristic scheme and also the CUSP scheme, provided that the coefficients of convective diffusion and pressure differences are correctly balanced. Scalar diffusion does not satisfy the first condition. In the case of the CUSP scheme (30) equation (36) reduces to
in the flux vector, which contain p H , and the state vector which contains pE. This may be remedied by introducing a modified state vector
Then one introduces the linearization fR
 f L = A h (WhR
.
 Wh'
z
Here A h may be calculated in the same wa as the standard Roe linearization. Introduce the weig ted averages defined by equation (27). Then Thus W R  W A is an eigenvector of the Roe matrix A R A , and $ is the corresponding eigenvalue. Since the eigenvaluesare U , U + c, and U  c, the only choice which leads to positive diffusion when U > 0 is U  c, yielding the relationship a * c = ( l + p >( C  U )
1
uH
^k
7
U
The eigenvalues of A h are U, A+ and A where
,o< U < c
Thus there is a one parameter family of schemes which support the ideal shock structure. The term p
Now both CUSP and characteristic schemes which preserve constant stagnation enthalpy in steady flow can be constructed from the modified Jacobian matrix A h [80]. These schemes also produce a discrete shock structure with one interiorpoint in steady flow. Then one arrives at four variations with this ropert ,which can conveniently be distin uished asthe E! and dCUSP schemes, and the E and Ifcharactenstic schemes. 4.5
Figure 5 : Diffusion Coefficients.
4.4.10
0
CUSP and Characteristic Schemes Admitting Constant Total Enthalpy in Steady Flow
In steady flow the stagnation enthalpy H is constant, corresponding to the fact that the energy and mass conservation equations are consistent when the constant factor H is removed from the energy equation. Discrete and semidiscrete schemes do not necessarily satisfy this property. In the case of a semidiscrete scheme expressed in viscosity form, equations (1 1) and (12),a solution with constant H is admitted if the viscosity for the energy equation reduces to the viscosity for the continuity equation with p replaced by p H . When the standard characteristic decomposition (28) is used, the viscous fluxes for p H which result from composition of the fluxes or and the characteristic variables do not have this roperty, and H is not constant in the discrete solution. fn practice there is an excursion of H in the discrete shock structure which represents a local heat source. In very high speed flows
P
Multidimensional Schemes
The simplest approach to the treatment of multidimensional problems on structured meshes is to appl the onedimensional construction separately in each mesh direction. On triangulated meshes in two or three dimensions the SLIP and USLIP constructions may also be implemented along the mesh edges [79]. A substantial body of current research is directed toward the implementation of truly multidimensional upwind schemes in which the upwind biasing is determined by properties of the flow rather than the mesh. A thorough review is given by Pailliere and Deconinck in reference [ 1321. Residual distribution schemes are an attractive a proach for triangulated meshes. In these the residual deined by the s ace derivatives is evaluated for each cell, and then distriiuted to the vertices with weights which depend on the direction of convection. For a scalar conservation law the weights can be chosen to maintain positivity with minimum cross diffusion in the direction normal to the flow. For the Euler equations the residual can be linearized by assuming that the parameter vector with components ,/&,@ui, and f i H varies linearly over the cell. Then
afj ( w )Aj&
axj
axj
where the Jacobian matrices A j = g are evaluated with Roe averaging of the values of w at the vertices. Waves
111
in the direction n can then be expressed in terms of the eigenvectors of njA3,and a positwe distribution scheme is used for waves in preferred directions. The best choice of these directions IS the subject of on oing research, but reliminary results indicate the possiiety of p+vinkkgh resolutlon of shocks and contact dmontlnuihes w cb are not aligned with mesh lines [132]. Hirsch and Van Ransbeeck adopt an alternative approach in whichthey directly conshuctduectionaldiffusive terms on structured meshes, with antidiffusion controlled by limiters based on comparisons of slopes in diffFent dIrections 601. They also show prouusing results in calculations o nozzles with multiply reflected oblique shocks.
c
4.5.1 Hi h Order Godunov Schemes, and Kinetic F l u
Sphing
A substantial body of current research is d e l e d toward the implementation of truly mul~dImensionalupwind schemes [59,135,101 Reference [132] provides a thorough review of recent velopments in this field. Some of the most impressive simulationsof time dependent flows with strong shock waves have been achieved with higher order Godunov schemes [1801. In these schemes the average value in each cell is updated by applying the integral conservation law using interfFe fluxes, predicted from the exact or approxlmate SOluhOn of a hemann problem between adiacent cells. A hieher order estimate of the solution IS then reconstructed h m the cell averages, and slope luniters are ap lied to the reconstruction. An example is the class ofessentially nonroscillatory (ENO) schemes, which can attain a very lugh order of accuracy at the cost of a substantial increase tn computational complexiy [32. 153, lS1,.152]. Methods b&d on reconstruction can also be unplemented on unstructured meshes 113, 121. Recently there has been an increasing interest in kinetlc flux splrtting schemes, which use solutions of the Boltzmann equatlcon or the BGK equation to predict the interface fluxes [42,36,45, 136, 1811.
For a cellcentered discretization (figure 6a)
2 is needed
at each face. The simplest procedure is to evaluate in each cell, and to average between the two cells on either side of a face [87]. %e resulting discretization does not have a com act stencil, and sup Its undamped oscillatory modes. a onediensionZalculation. for U + 2U.tUI2 example, would be discretized as * ' 4kz . In order to produce a compact stencil may be estimated From acontrol volume centered on each face,using formulas (38) or (39) 11441. This is com utationally expensive because the number of faces is mu& larger than the number of cells. In a hexahedral mesh with ti large number of vermes the number of faces approaches three times the number of cells. This motivates the [email protected] of dual meshes for the evaluatlon of the velocity denvatlves and the flux balance as sketched in figure 6. The figure shows both
&
2
R
2
dk
6a: Cellcentered scheme. uij evaluated 6b: Cellvertex scheme. at vertices ofthe primaryoij evaluated at cell centers of the primary mesh mesh
Fi ure 6 Viscous discretizations for cellcentered and ceivertix algorithms. Tbe discretization of the viscous terms of the Navier Stokes equations requires an approximation to the vclocity derivatives in order to calculate the tensor u,j, defined by uation (3). Then the viscous terms may be included i n x e flux balance (4 In order to evaluate the derivatives one may apply the auss formula to a control volume V with the boundary S
2
where nj is the outward normal. For a teuahedral or hexahedral cell this gives
cellcentered and cellvertex schemes. The dual mesh connects cell centers of the orimarv a~. ,mesh. If . .there is ~. kink in the ~~&a&mesh:&ule'du&xlls should be formed by assemdng contiguous fractlons of the neighboring primary cells. On smooth meshes comparable results are obtained b either of these formulations [114,115, 1071. If the mest has a kink the cellvertex scheme has the advantage that the derivatives are calculated in the interior of a regular cell, with no loss of accuracy. A desirable property is that a linearly varying velocity distribution, as in a Couette flow, should produce a constant stress and hence an exact stress balance. This roperty is not necessarily satisfied in general by finite digerence or finite volume schemes on curvilinear meshes. The characterization kexact has been proposed for schemes that are exact for polynomials of degree k. The cellvertex finite volume scheme is linearly exact if the derivatives are evaluated by equation (39). since then is exactly evaluated as a constant, leading to constant viscous stresses uij, and an exact viscous stress balance. This remains true when there is a kink in the mesh, because the summation of constant stresses over the faces of the kinked control volume sketched in figure 6 still yields a perfect bowbalance. The use of equation (39) to evaluate ever, requires the additional calculation or storage of the in each cell, whereas equation nine metric quantities (38) can be evaluated from the same face areas that are used for the flux balance. In the case of an unshuctured mesh, the weak form (6) leads to a natural discretization with linear elements, in ~~~~~~~
2
2
where a, is an estimate of the avera e of U, over the face. If u varies linearly over a tetmiedral cell this is exact. Alternatively, assuming a local transformation to computational coordinates tj one may apply the chain
2,
e
can be evaluated Here the transformation derivatives by the same finite difference formulas as the velocity In this case is exact if U is a linearly derivatives varying functlon.
e
e
2
112
which the piecewise linear approximation ields a constant stress in each cell. This method yieldys a representation which is globally correct when averaged over the cells, a result that can be proved by energ esumates for elliptic problems [ 1641. It should be. n o t d however, mat it yields formulas that are not necessarily locally consistent with the differential equations, if Ta lor seriesexpansions are substltuted for the solution at tie vertices [email protected] in the local stencil. Figure 7 illustrates the mscretlzatlon of the Laplacian uZz+ uyywhich is obtained with linear elements. It shows a particular triangulation such that the approximation is locally consistent with uzz + 3uyy. Thus the use ofan irregular uian lation in the boundary layer may significantly degrade tre accuracy.
h
b
h
Figure 7 Example of discretization..u + uyuon a triangular mesh. The discretizationis locally equivalent to the approximation U,,=*, 3~,,="d6hU,.+~~ .
Time Stepping Schemes If the space discretization rocedwe is implemented separately, it leads to a set ofcoup!ed ordinary differentlal equatlons, which can be written in the form
4.7
dw dt
 + R(w)=O, where w is the vector of the Bow variables at the mesh points, and R(w) is the vector of the residuals, consisting of the flux balances defined by the space discretization scheme, together with the added dissi ative terms. If the obiective is simolv to reach the stedv state and details
derivatives are calculated from known values of theflow vaiables,at the, beginning of the time step, or an kplicit scheme, UI whch the formulas for the s ace denvatlves include as yet unknown values of the &w variables at the end of the time ste , leading to the need to solve coupled equations for &e new values, p e pemsstble time step for an explicit scheme IS h t e d by the CourantFnednchsLewy (a) condition, which states that a differencescheme cannot be a conver ent and stable a prodation unless its domain of depenfence contains $e d o m m of dependence of the corresponding Merential eouation. One can anticioate that unolicit r~~~~~ schemes will geld convergence in a sm'aller number of time ste s because the timestep is no longer constrainedby the Ck condjtion. Imphcit schemes will be. efficient, however, only if the decrease in the number of ume steps outweighs thehcrease in the corn utational effort per time step consequent u p n the n d t o solve coupled equations. The prototype lmplicit scheme can be formulated by estimating at t + pAt as a linear combination of R (w") and ~~
~
~
~~~~~~~~~~~~
~
~
~~~~
~
R(w"+'). The resulting equation
can be linearized as
If one sets p=l and lets At + m this reduces to the Newton iteration which hm been successfully used in twodimensional calculations 1173, 50 In the threedimensional case with, sa an N x ". x N mesh, the bandwidth of the matrix tkat must be inverted is of order N'. Direct inversion requires a number of operations proportional to the number of unknowns multlplied by the s uare of the bandwidth of the order of N7.This is rohhtive, and forces recourse to either,an approximate Factorization method or an iteranve solutlon method. Alternatin direction meth+, which introduce factors corresponiing to each coordinate, are widely used for structured meshes [17, 1371. They cannot be implemented on unstructured tetrahedral meshes that do not contain identifiable mesh dxections, although other decompositions are possible [log]. If one chooses to adopt the iterative solutlon techruque, the rinci al alternatives are variants of the GaussSsidel a n i Jacogi methods. A symmetric GaussSeidel metpod with one iteration per tlme step is essentially eqtyalent to an approximate lowerupper (LU)factorizatlon of the imphcit scheme [86,125,31,184]. On the other hand, the Jacobi method with a fixed number of iterations per time step reduces to a multista e explicit scheme, belongin to the en eral class of fun e Kutta schemes [33 Sciemes o?thiI type have rovdv& effective for wide varie of prob lems, and %ey have the advanta e that they c a n L applied equally easily on both structurdand unstructuredmeshes [84,67,69, 1451. If one reduces the linear model problem corresponding to (40) to an ordinary differential equation by substitutinga Fourier mode fi=e'PzJ, the resultin Fourier symbol has t proportional to fhe wave speed. and an imapnary a a negatlve reaf'part roportional to the diffusion. n u s the tune stepping sc?I erne should have a stability region which contams substantial intervals of both the negative real axis and the hagin? yis.To achieve this it pays to treat the convectlve an dssipatlve terms in a distlnct fashion. Thus the residual is split as
.
R(w)=Q(d + D(ur), where Q(,,,)is the convective and D(w) the dissiative art. Denote the time level nAt by a superscri t n h e n g e multistage time stepping scheme is f o r m u f k as ,,,(n+i,o)
=
,,n
113 chosen to increase the stability interval along the negatlve real axis. These Schemes do not fall the mework of RungeKutta schemes, and they have much larger stability regions [69]. live schemes which have been found to be particularl effective are tabulated below. The first is a fourstage scieme with two evaluations of dissipation. Its coefficients are
has to be transferred back to grid k  1 with the aid of
..
an intemolationoveratorL .7 . .,k .. With uroverlv  . ovtimized coefficikntsmultiitagetimestepping schemes can be very 'd P'JCess. A WcYcle of effiC*ent,fiven of the the type illustrated In [email protected](lfProves to be a Particularly
The second is a fivestage scheme with three evaluations of dissipation. Its coefficients are
4.8 Multigrid Methods
4.8.1 Accelerntion of Steady Flow Calculations
Radical improvements in the rate of convergence to a steady state can be realized by the multigrid timestep ing t e c h ue The concept of acceleration by the intrduction o? mhtiple grids was first proposed b Fedorenko [48]. There is by now a fairly welldevegped theory of multi 'd methods for elliptic e q d o n s based on the concept g a t the updating scheme acts as a smoothing o p erator on each grid [24,53]. This theory does not hold for hy rbolic s stems. Nevertheless, it Seems that it ou ht to possibz to accelerate the evolution of a hyperb3ic system to a steady state by usin large [email protected] s on coarse ds so that disturbances wilfbe more rap&' expelled g o u g h the outer boundary. Vanous multl$d timestepping schemes designed to take advantage o t h ~ effect s have been . proposed [124,65,55,71,29.6,57,83,93]. . One can devise a multi 'd scheme using a sequence of independently g e n e r a t T c o a y meshes, by ehminating alternate pomts in each coordinate %cpon. In order to give a precise descnptlon of the mule d scheme, subscripts may be used to indxate the FSeveral transfer operationsneed to be defined. First solution vector on grid k must be initialized as
E
%
8c:5 Levels. Figure 8: Multigrid WcYcle for managing the grid Galculation. E, evaluate the change in the Bow for one step; T,transfer the data without updating the solution. effectivestrate formanagingthe worksplitbetween the meshes. In a g w d i e n s i o n d Case the number of cells is reduced by a factor of eight on each coarser grid. On examination of the figure,it can thereforebe seen that the work measured in u t s correspondmg to a step on the fine grid is of the order of 1 + 218 + 4/64 +
... < 413,
and consequently the very lar e effectivetime step of the complete cycle costs only slig tly more than a single time step in the fine grid.
%
4.8.2 Multigrid Implicit Schemes for Unsteady Flow
where wk1 is the current value on grid k  I, and !fk,kI is a transfer operator. Next it is necessary to transfer a residual forcing function such that the solution grid k is driven hy the residuals calculated on grid k  1. This can be accomplished by setting Pk=Qk,kiRki
(wkI)
 R I [Wp)] ,
where Qk,kI is another transfer operator. Then Rk ( W k ) is replaced by Rk (Wk + P k in the time step ing scheme. Thus. the multistage scheme is reformulateaas
w!')
=
Pk]
Dtw""
+ R(w"+') =O.
Here Dt is a kth order accurate backward difference operator of the form
...
... W k('I')
tup)  alAttn: [Rp)+
T h e dependent calculations are needed for a number of important ap hcatlons. such as Buner analysis, or the analysis of the low pas1 a helicopter rotor, in which the stability linut ofan explicit scheme forces the use of much smaller time steps than would be needed for an accurate simulation. In this situation a multigrid explicit scheme can be used in an inner iteration to solve the equations of a fully implicit time stepping scheme [74]. Suppose that (40) is approximated as
=
W p )  cY,+l&
[ , : I '
+ Pk]
.
Tbe result wi'") then provides the initial data for grid k + 1. Finally, the accumulated correction on grid k
where
Aw"+l
w "+I w".
114
Applied to the h e a r differential equation
dw dt
=(2W
the schemes with k = l , 2 are stable for all aAt in the left half plane (Astable). Dahlquist has shown that Astable linear multistep schemes are at best second order accurate 381. Gear however, has shown that the schemes with < 6 are stim stable [49], and one of the higher order sckmes may o er a better compromse between accuracy and stability, depending on the application. Equation (40) is now treated as a modified steady state problem to be solved by a multigrid scheme using variable local time steps in a fictitious time to. For example, in the case k=2 one solves
I
4
aw
=R*
at*
(w)
in the different coordinate directions. The need to resolve the boundary layer generally forces the intiduction of mesh cells with very high aspect ratios near the boundand these can lead to a severe reduction in the rate yionvergence to a steady state. Pierce has recently obtained impressiveresults using diagonal and blockJacobi preconditloners which include the mesh intervals [ 1331. An alternative approach has recently been proposed by Ta'asan 1681, in which the equations are wntten in a canonical form which se arates the equations describing acoustic waves from ose describing convection. In terms of the velocity components U,v and the vorticity w , temperature T,entrop s and total enthal y H, the eauations describinn s t e d twodimensional &w can be
J:
[
,
where where and the last two terms are mated as k e d source terms. The first term shifts the Fourier symbol of the equivalent model problem io the left in the complex plane. While this pr6motes stability, it may also reuirea limit to be imposed on the magnitude of the local time step At* relative to that of the implicit time step At. This may be relieved by a oint im licit modification of the multistage scheme [ fl9].In %e case of problems with moving boundaries the e uations must be modified to allow for movement and delormation of the mesh. This method has proved effective for the calculation of unsteady flows,that mi ht be associated with,wing flutter and also in the ciculation of unsteady incompress!&?lows 1181. It has the advantage that it can be added as an optioh tdacom uterprogrdwhich uses an ex licit multi 'd scheme, alfowing it to be used for the e d i e n t calcuEon of both steady and unsteady flows. 4.9 Preconditioning
Another way to im rove the rate of convergence to a steady state is to muPtiply the space derivatives in equation (1) by a reconditioning matrix P which is designed to equalize $e eigenvalues, so that a~ the waves can advanced with optlmal tune steps. A symmetnc preconmtioner which ualizes the eigenvalues has been proposed by Van Leer3021. When the equations are written in streamaligned coordinatesthis has the form
r &MZ
 IPM o o o &+l 0 0 0
1
where
Turkel has pro osed an asymmetric preconditioner which has also provdeffective.particularly for flow at low Mach numbers [ 1721. The use,of these preconditionencan lead to instability at stagnauon points where there is a zero eigenvalue which cannot be equalized with the eigenvalues tc. The preconditionem of Van Leer and Turkel do not take account of the effect of differences in the mesh intervals
D3
Q
=
a az
UU
a
= utv
az
a ay 8
ay
Here the first two quatio,ns describe ?,elliptic system if the flow is subsomc, wbde the remaning equations are convective. Now se arately optimized mulugrid procedures are used to the two sets of equations, which are essentially decoupled.
sage
4.10 High Order Schemes and Mesh Rehement The need both to improve the accuracy of computational simulations and to assure known levels of accuracy is the focus of ongoing research. The main routes to improvin the accuracy are to increase the order of the discrete scfeme and to reduce the mesh interval. Hi h order difference methods are most easily implemente3 on Cartesian, or at least extreme1 smooth gnds. The expansion of the stencil as the ordkr is increased leads to the need for complex boun conditions. Com act schemes keep 281. On simple the stencil as s a a s p s i b l e [ 140, domains. spectral me ods are articularly effective, especially in the case of periodic Lundary conditions, and can be used to produce exponentially fast convergence of the mor as the mesh interval is decreased 1127,271. A compromise is to divide the field mto subdomains and introduce high order elements. This approach is used in the spectral element method [92]. High order difference schemes and s ctral methods have proven articular~yuseful in direct Xvierstokes simulauons o?transient @d turbulent flows. High order methods are also beneficial in computational aeroacousucs. where it is desired to track wavds over long distanceswith minimum error. If the flow contains shock waves or contact discontinuities, the EN0 method may be used to construct high order nonoscillatory schemes. In multidmensional flow simulations, global reduction of the mesh interval can ,be prohibitively expensive, motlvatine. the use of adaouve mesh refinement urocedures whichduce the local mesh width h if there is' an indication that the error is too large [21,39, 109,61,138. 1031. In such hrefinement metfiods, simple error indicators
lb,
115
such as local solution gradients may be used. Alternatively, the discretization error may be estimated b comparing quantities calculated with two mesh wid& say on the current mesh and a coarser mesh with double the mesh interval. Procedures of this kind may also be used to provide a posteriori estimates of the error once the calculationis completed. This kind of local ada tive control can also be applied to the local order of a h e element method to produce a prefinement method, where p represents the order of the polynomial basis functions: Finally, both h and p refinement can be combined to produce an hp method in which h and D are locallv ootimized to vield a solution with minimum ermr [ 126].rSuch m e t h h can achieve exponentially fast convergence, and are well established in computational solid mechanics. ~~~~~
~
. ~ ~  ~
~~~~~
~~~
5. CURRENT
LATION
STATUS OF NUMERICAL SIMU
This section presents some representative numerical results which con& the pro ernes of the al orithmswhich have been reviewed in theyast section. &ese have been drawn from the work of the author and his associates. The also illustrate the kind of calculation which can be pert?&ed in an industrial environment, where rapid turn around is mportant to allow the quick assessment of design changes, and computational costs must be Limited. 5.1 Onedimensional shnck
In order to ve
the discrete structure of stationary shocks, calcu atlons were erformed for a onedimensional problem with initial &ta containing left and right states compatible with the Rankine Hugniit con$tions. An intermediate state consisting of e anthmetlc average of the left and ri ht states was introduced at a single cell in the center offhe domain. With this intermediate state the svstem is not in eouilibrium. and the time de ndent equatlons were solved to find an equilibrium s o h o n with a stationary shock wave separatlng the left and right states. Table 1 shows the result for a shock wave at Mach 20. This calculation used the HCUSP scheme, which allows a solution with constant sta na tion enthal y. with the l i t e r defined by equation &3), and q=3. &e formulation is described in detad in reference [SO]. The table shows the values of H, p , M and the entropy s=log  log A perfect one point shock structure is displayed. The entrop is zero to 4 decimal places upstream of the shock, e&bits a slight excursion at the mterior oint, and is constant to 4 decimal places downstream ofthe shock. It may be noted that the mass, momentum and energy of the initial data are not compatible with the final uilibrium state. According to conservatlon arguments e total mass, momentum and energy must remain constant if the outflow flux fR remains ual to the inflow flux fL. Therefore fR must be al1ow3 to vary according to an appropriate outflow boundary condition to allow the total mass, momentum and energy to be adjusted to values compatible with equilibrium.
9.
~~~~~~
~~~~~~
~~~~
~~
A
~~~~
~~~
~
~~~~~
~~~~
~~~~~
5.2 Euler cnlculationsfor Airfoils and Wings
The results of transonic flow calculations for two well known airfoils, the RAE 2822 and the NACA 0012, are presented in figures (2225). The HCUSP scheme was again used. The Limiter defined by equation (23)was used with 4=3.The 5 stagetime steppin scheme (42)was augmented by the mulhgrid scheme &scribed in section, 4.2 to accelerate convergence to a stead state The equatlons were discretized on meshes with 6topoiogy extending out to a radius of about 100 chords. In each case the calculations were performed on a se uence of successively finer meshes from 40x8 to 320x84 cells, while the multlgrid cycles on each of these meshes descended to a coarsest mesh of 10x2 cells. Fi 22 shows the inner parts of the 160x32 meshes for %wo airfoils. Figures 2325 show the h a l results on 320x64 meshes for the RAE 2822 airfoil at Mach .75and '3 angle of attack, and for the NACA 0012 airfoil at Mach .8 and 1.25" angle of attack, and also at Mach .85 and 1' angle of attack. In the pressure distributions the pressure coefficient C P = m is dotted with the neeative (suction) uoward. so ,~~~ ,oressures r thh the u p p c r c y e g resentsthe flow overthe tipper side of a lifting aufoil. &e convergence lustones show the mean rate of change of the density, and also the total number of supersonic points in the flow field, which provides a useful measure of the global conver ence of tiansonic flow calculations such as these. In ea& case the convery e history is shown for 100 cycles, while the pressure smbuGon IS displayed after a suflicient numkr of c cles for its convergence. The pressure dtstnbutlon of& RAE 2822 airfoil conver ed in only 25 c cles. Convergence was slower for thekACA 0012 aidil. In the case of flow at Mach .8 and 1.25' angle of attack, additional c cles were needed to damp out a wave downstream of weak shock wave on the lower surface. As a further check on accuracy the dra coefficient should be zero in subsonic flow. or in shock f r e transonic ~ flow. Table 2 shows the corn' uted drag coefficient on a sequence of three meshes three exam les The first two are subsonic flows over the RAE 282fand NACA 0012 airfoils at Mach .5 and 3" angle of attack. The third is the flow over the shock free Kom airfoil at its design point of Mach .75and 0' angle of attack. In all three cases the dra coefficient is calculated to be zero to four digits on a l6&32 mesh. ~~~~
~~~~~~~
~
de
~~~
~~~~~~~~~~~
P,
(3).
Table 2: Drag Coefficient on a sequence of meshes
%
n
283.5ooO
283.5000 283.4960 283.4960 283.4960 283.4960
NI
S
zmRxT(I
0.0000
1 .0000 20.0000 0.0000 1.0000 20.0000 0.0000 307.4467 0.7229 40.3353 466.4889 0.3804 37.6355 466.4889 0.3804 37.6355 466.4889 0.3804 37.6355
Table 1: Shock Wave at Mach 20
and 3.06' angle Gf attack. This again verifies the nonoscillatory character of the solution, and the sharp resolution of shock waves. In this case 50 c cles were sufficient for convergence of the pressure dismiutions. Figure9showsacalculationoftheNorthropYF23byR.J. Busch. Jr., who used the author's K O 5 7 code to solve thc Euler equations 1261. Although an inviscid model of the flow was used, it can be seen that the simulationsare in ood a ement with wind tunnel measurements both at hach with angles of attack of 0.8 and 16 degrees, and at Mach 1.5 with an le, of attack of 0, 4 and 8 deees. At a hgb an le ofakick the flow separates from g e leading edge, anithis example showsthat in situations where the point of separation is fixed, an inviscid model may still produce a useful predictlon. Thus valuable in
.6
116
I No. ot Nodes 11
formation for the aerodynamic design could be obtained with a relatively inexpensive computational model.
SecondslCycle [ Speedup
I
Table 3: AIRPLANE Parallel Performance on the SE, MD11 Model
~
...
~
~~~
~~~~~~~~
~
~
~~~~
~~~~
throu b the engine nacelles, using 3%8407mesh points of
2 1 d 6 6 tetrahedra. This calculation takes 4 hours on an IBM 590 workstation. A parallel version of the code has been developin collabohion with W.S. Cheng, and !he same calcu anon can be performed in 20 m u t e s using 16 mcessms of an IFJM SPZ. The parallel speedup for the%iDll is shown in table 3.
I

Figure 9 Comparison of Ex erimental and Com uted Drag Rise Curve for the YF2!3 (Supplied by R. J.\usb Jr.)
Figure 11: Pressure Contours and Sonic Boom on a Representative HSCT Configuration
Vioeous Flow Calculations The next figures show viscous simulations based on the solution of the Reynolds avcra ed Navier Stokes equations with turbulence models. I ure 13 shows a twodimensional calculation for the 2822 airfoil by L. Martinelli. The vertical axis re sents the ne ative res sure coefficient, and there is a sKck wave halfway $ong the upper surface. This exam le confirms that in the absence of significant shock in uced separation, simulations performed on a sufsciently fine mesh (with 5 12 x 64 cells) can produce excellent agreement with ex rimental data. Fi e 21 shows a simulation of the ZDonneUD o u g l a s ~ 8performed by R.M.Cummin s, YM Rizk, L.B. Schiff and N.M. Chaderjian at NASS A& [37 They used a multiblock mesh with about 900000 mes points. While this is probably not enough for an accurate. quantitative prediction,.the agrFment with both the expenmental data and the vlsualizatmn are quite good. Figure 14 shows an unstead flow calculation for a itchin airfoil performed b Alonso usin the code which he 'ointly d v e l y d wi$ L.%artinelli and the author [41. uses them tlmd un~licitscheme described in Section 3.7.2 which alhws thi number of time steps to be reduced from several thousand to 36 per pitching c cle The agreement with experimentaldata is qLIitego0J . 5.3
& '
B
k
Figure 1 0 Com arison of Experimental and Calculated Results for a H S f f Configuration The next fi ures show the results of calculations using the AIRPLAN% code developed by T.J. Baker and the author, to solve theEuler uations onan unstrucnurdmesh. This proyides the flexit& to treat arbitrarily complex configuratlons without the need to spend months developing an
h0h
As
f.
117
4
Figure 12: Computed Pressure Field for a McDonnell Douglas MDl 1 5.4
Ship Wave Reaistauce caIculations
Figure; 1517 show the results of an application of the same multigrid h i t e volume techniques to the calculation of the flow past a naval fri~te,,usinga code which was developed by J. Farmer, L. artrnelh and the author [47]. The mesh was adjusted during the course of the calculation to conform to the free surface in order to satisfy the exact nonlinear boundary condition, while ,?rtificial compressibility was used to treat the incompressible flow equations. 6. AERODYNAMIC SHAPE OpTIMlzATION
Figure 13: noDimensional 'hbulent Viscous Calculation (by Luigi Martinelli)
6.1 Optimization and Design
Traditionall the rocess of selecting design variations has been carridout I! & ! d and error, relying on the intuition and experience orthe desi ner With currently available equipment the turn mount foi numerical simulations is becomin so rapid that it is feasible to examine an extremely &ge number of variations. It is not at likely that re ated trials in an interactive desi n and analysis proce& can lead to a truly optimum &sign. In order to take full advanta e of the possibili of examinin a large design s ace e numencal simu auons need to%e combined w i g automatic search and optimization procedyes. This cap lead to autopatic design meth.?ds which will full realize the potenual unprovements in aerodynamic ehciency. The simplest ap roach to optimization is to define the eometry througg a set of design payneters, which may, for example, be the weights ai applied to a set of shape functions b, (z) so that the shape is represented as
7.
ti
Then a cost function I is selected which might, for example, be the drag coefficient or the lift to drag ratio, and I is regarded as a function of the parameters CY,. The sensitivities may now be estimated by making a small variation CY. in each design parameter in turn and recalculating the flow to obtain the change in I . Then
%
The gradient vector may now be used to determine a direction of im rovcient.  m e simplest procedure is to make a step ine! negative gradient direction by setting CY"+'=CY"
 A60,
so that to first order
I + 6I=I  &=I ap aa
ar aCYaCY
 A.a r T
More so histicated search rocedures may &used such as quasidwtonmethods, attempt to estunate the second derivative of the cost function fromchanges in the Merit in successive optimization steps. These methods also generally introduce line searches to find the minimum in the search direction which is defined at ~. each step. The main disadvantage ofthis approach is the need for a number of flow calculations proponional to the number of design variables to estimate the gradient. The computational costs can thus become prolubitive as the number of design variables is increased. An alternative approach is to cast the desi n problem as a search for the shape that will generate the &red ressure distribution. This approach recognizes that the iesigner usually has an idea of the the kind of pressure distnbution that will lead to the desired performance. Thus. it is useful to consider the inverse problem of calculating the sha that will lead to a given pressure distribution. The megod bas the advanta e that only one flow solution is required to obtain the &sired design. Unfortunately, a bysically realizable shap may not necessarily exist,,unPess the pressure hstnbutlon satlsfies certllln consmmts. Thus the problem must be very carefully formulated,otherwise it may be ill posed. The difficulty that the target pressure may be unattainable may be circumvented by treating the inverse problem as a special case of the optimization problem, with a cost function which measures the error m the solution of the inverse problem. For example, if pd is the desired surface pressure, one may take the cost function to be an integral over the the body surface of the square of the pressure
&
with
~~~~~~
~~
~~~~~~~~
118
Figure 1 6 Contours of Surface Wave Elevation Near the Transom Stem
Figure 14 Mach Number Contours. Pitching Airfoil Case. Re=l.O x IO6,M,=0.796,K,=0.202. error,
or possibly a more general Sobolev no? of the pressure error. This has the advantage of converung a possibly ill posed problem into a well posed one. It has the disadvaniage h i t it incurs the computational costs associated with optimization procedures. 6.2 Application of Cnntml Theory
In order to reduce the computational costs, it turns out that there are advantages in formulating boththe inverseproblem and more eneral aerodynmc problems w i h n the framework of &e mathematlcal theorv for the control of s stems overned by partial differenual e uations [lo5 lwing, for exayle. is a,deviceto producegift by control: line the flow, an its design can be regarded as a problem in &e optimal control ofthe flow equations by variation
Figure 17: Pressure Contours in the Bow Region of the shape of the boundary. If the boundary shape is regarded as arbitr within some requirements of smoothness, then the fuygenerality of shapes cannot be defined with a finite number of arameters, and one must use the concept of the Frechet aerivative of the cost with respect tn a function. cannot~be deter.... Clearlv.  such, a derivative , ~ ~ mined directly by filute differences of the design arameters because there are now an infinite number orthese. Using techniques of conml theory, however, the,gradient can be determined indirectly b solvin an adjoint e ua tion which has coefficients defiYned by he soluuon 09th; flow equations. The cost of solving ihe adjoint equation is comparable to that of solvin the flow e uations. Thus the gradient can be determindwith roughy the computational costs of two flow solutions, indegqdently of the number of design variables, which may infinite if the boundary is regarded as a free surface. For flow about an airfoil or wing, the aerodynamic pro erties which define the cost function are functions of flowfield variables ( w ) and the physical location of the boundary, which may be represented by the function 3, say. Then ~
~~~
~~~~~~
~~~~
~
~~~~~~~~~~
~~~~
I=I(w,rn, and a change in 3results in a change ____
Figure 15: Contours of Surface Wave Elevation for a Combatant Ship
IT
61=6w aw
aiT
+ 63, 83
(43)
in the cost function. Using control theory, the governing
~
119
equations of the flowfield arc introduced as a consuaint in such a way that the final ex ression for the gradient does not requue reevaluation oI%e flowfield. In order lo achieve this 6w must be elimina!ed from (43). Su pose that the govenungequauon R which expresses the J p e n dence of w and 7 within the flowfield domain D can be written as R (w,7 )=O. (44) Then 6w is determined from the equation (45) Next, introducing a Lagrange Multiplier $, we have 6I
=
aIT
arT
aw
a3
+
6w+~5$~
=
$:[
T
[8;;;])6w+[$$ aR
T
[])63. aR
aF
Choosing (I, to satisfy the adjoint equation
the first term is eliminated. and we find that
61=867,
(47)
where
The advantage is that (47) is independent of 6w, with the result that the gradient of I with respect to an arbitraty number of desien variables can be determined without the need for additi&fiowfield evaluations. In the case that (44)is a artial differential equation, @eadjoint equation (46) is a partial Merenual equauon and appropnate boundary conditions must be deteimined. After makine a s t e ~in the neeative eradient direction. the gradient Fan be iecalculadand thi process repeated to follow a path of steepest descent unul a minimum is reached. In order to avoid violaung consmamts, such as a minimum acceptable wing thickness, the gradient ma be projected into the allowable subspace within whic; the constraints are satisfied. In this way one can devise procedures which must necessarily converge at least to a ocal minimum. and which can be accelerated bv the use 7  ~ of more sophisticated descent methods such as conjugate %radmdor quasiNewton alfri%s. There is the possiility of more than one loc mmimum. but in any case the method will lead to an improvement ovei the original design. Furthermore, unlike the traditional mverse algor i h i s , any measure of performance can be used as The cost funchon. In reference [72] the author derived the adjoint equations for transonic flows modelled by both the potenhal flow equation and the Euler equations. The theory was developed in terms of artial differential equations, leading to an adjoint partiJ differential equation. In order Jo obtain numerical solutions both the flow and the adjomt equations must be discretized. The control theory might be a plied directly to the discrete flow uations whch from the numerical amroximation3the flow eauations hy finite element, fin& volume or finite differ&cc procedures. This leads directly to a set of discrete ad'oint equations with amatrix which is the trans ose of the hcobran mahx of the full set of discrete noninear flow q u a tions. Onathreedimensional mesh with indicesi, ,kthe individual adjoint equationsmay be denved by cokxting together all the terms multiplied by the variation 6 W , j , k of the discrete flow variable W , , j . k . The mulling discrete adjoint equationsrepresent a possible discretization of the ~~~~~~~~~
~~
~~~~
go
~~~~~
~~~~
~
~~~~~
~~~~
~~
d
red
~
~
~
~~
~~~~
~~~
ad'ointpartial differentialequation.Ifthese uations are solved exactly they can provide an exact grzent of the inexact cost function whch results from the discretization of the flow equations. On the other hand any consistent discretization of the adjoint partial differenual equation will yield the exact adient in the limit as,the mesh is refined. The tradeohetween the compleuty of the adjoint discretization. the accuracy of the resulting estimate of the gradient, and its impact on the computauonal cost to approach an optimum solution is a subject of ongoing research. The me optimum shape belongs to an infinitely dimensional space of design arameters. One motivation for developing the theory !or the partial differential q u a tions of the flow is to provide an indication in principle of how such a solution could be ap roached if sufficient computational resources were ava&ble. Another motiyation is that it highlights the possibilit of generatin dl posed formulauons of the problem. &or example, i one attempts to calculate the sensitivity of the pressure at a particular location to changes in the boundary shape, there is the ossibility that a sha e modification could cause a shocf wave to pass over &t location. Then the sensitivity could become unbounded. The movement of the shock, however, is continuous as the shape changes. Therefore a uantity such as the drag coefficient, which is detennine8by integrating the pressure over the surface, also depends continuously on the shape. The adjoint equation allows the sensihvity of the drag coefficient to be determined without the ex licit evaluauon of pressure sensitivities which would be &posed. The discrete adjoint equations, whether they are derived directly or by discretization of the adjoint partial differential equation, are linear. Therefore they could be solved by direct numericalinversion. The cost of direct inversion can become prohibihve, however, as the mesh is relined, and it becomes more efficient to use iterative solution methods. Moreover, because of the similarity of the adjoint y t i o n s to the flow equations, the s F e iterative metho s wluch have been proved to be efficient for the solution of the flow equations are efficient for the solution of the adjoint equations. The control theory formulation for optimal aerodynamic desi n has roved effective in a vanety of ap lications [73.57, 144. The adioint equations have also &en used by Ta'asan, Xuruvilaand Sdas [167 who have implemented a one shot approach in which e constraintrepresented bv the flow euuations ~ is onlv muired to be ~ satisfied by the fihal converged solution, and computational costs are also reduced by applying multi 'd techniques to the geometry modifications as well as &?solution of the flow ind ad'oint equations. pironneau has studied the use of controi theory for optimal sha desi n of systems governed b elliptic equations [ 1 6 , antmore recently the Navierstokes equations, and also wave reflection problems. Adjoint methods have also been used by Baysal andEleshaky [16].
j,

6.3 ThreeDimensionalDesign using the Euler QuatiOnS
In order to illustrate the application of control theory to aerodynamic design roblems, this section treats the case of thrwdimensionafwing design usin the inviscid Euler equations as the mathematical mode for compressible flow. A,transformation to a bodyfitted coorqinatesystem will be introduced, so that vanauons in the wing shape mduce corresponding variations in the computationalmesh. Thus the flow is determined by the soluhon of the transformed equation (5). Let
!
and
Q=JK~.
120
The elements of Q are,the coefficients of K,and in a finite volume discretizahongey are just the face areas of the computational cells projected in the 2 12 2~. and 2 3 directions. Also introduce scaled contravanant velocity components ui=Qijuj. The msformed equations can now be written as
aw + =o OF+ at
where
R8
The weak form of the equation for 6w in the steady state becomes
where
~F,=C;~W + 6Qijfj,
which should hold for any differential test function 4. This equation may he added to.the variation in the cost function, wluch may now be wntten as
W=Jw
and
Assume now that the new computational coordinate system conforms to the wing in such a way that the wmg surface BW is represented by E 0. Then the Bow is determined as the steady state solition of equation (48) subject to the flow tangency condmon
(53) On tpe win surface Bw,nl=n,=O and it follows from equahon (48) that
U 2 4 onBw. (49) At the far field boundary Bp, conditionsare specified for incoming waves. as in the twodimensional case, while outgoing waves are determined by the soluhon. The weak form of the Euler equationsfor steady Bow can be written as
(54)
Since the weak equation for bw should h r an arbitrary choice of thitest vector 4, we are free to choose 4 to
where the test vector 4 is an arbitrary differentiable function and n. is the outward normal at the boundary. If a differentiahe solution w is obtained to this equation, it can be integrated by parts to gwe
and since this is m e for any 0, the differential form can be recovered. If the solution is discontinuous, equation (5%may be integrated by parts se arately on either side of e discontinlutyto recover the sock jump condtions. Suppose now that it is desired to control the surfacepressure h varying the win shape. It is convenient to retam x Jcomputational omm. Vafiatons in the shape a6 then result rn correspondmg vanahons in the mapping derivatives defined by K. Introduce the cost function
f
simplify the resulting expressions. Therefore we set 4=$, whe? the costate vector 11, is the solution of the adjoint equahon (55)
At the outer boundary mcoming characteristicsfor $ correspond to outgoing characteristics for 6w. Consequently one can choose boundary conditionsfor 11, such that n.11,TC,6w=0. Then if the coordinate transformation is such that 6 negligible in the far field, the only remaining boun tern 1s
/Lw
% IS
11,T6Fz @[email protected]
Thus by letting $ satisfy the boundary condition, where p d is the desired pressure. The design problem is now treated as a control oblem where the control function is the wing shape, wgch is to be chosen to minimize I subject to the constraintsdefined by the flow equations (4850). A variation in the shape will cause a variation 6p in the pressure and consequently a variation in the cost function
Qz& + Q& + [email protected][email protected]  P d ) on Bw, we find finally that bI

=
(56)
J, %
6Qijfjd’D
//ca~,,$~
+ 6Qn$3 + Q d ~ 4 ) ~ C y i d h (57) .
BW
Since p depends on w through the equation of state (2), the variation 6p c,an be determined from the variation 6w. Define the Jacobian mamces
A convenient way to treat a wing is to introduce sheared parabolic coordinates as shown in figure 18 through the transformauon z =
1
20
(C)+ p (C) {E2  (II+ s(E, C 9 )
Y
= ~o(C)+a(C)E(II+S(E,C))
2
=
C.
121
18b: (,qPlane.
18a: 5 , yPlane.
Figure 18: Sheared Parabolic Mapping. Herez=zl y=zz. z = q are the Cartesian coordinates, and E andn + 3 correspond to parabolic coordinates generated 6y the'mapping 2
1 + iy=zo + iyo + a(()
2
{t + i (q +
at a fixed span station C. 50(0,and yo (C) ?e, the,coordinates of a sineular line which IS swem to he lust inside the leadin of a swept win , wihe a (C jis a scale factor to d o w for spanwise variations.
$,
chad
We now treat S C) $e control. Substitution of these formulas yields e vanatlon in the form
Inde ndent movement of the boundary mesh points coulrproduce discontinuities in the designed shape. In order to r y e n t ti$s @e,gradientmay be also smoothed. Both exp cit and imphcit smoothing procedures are useful. Sup sethatthemovement ofthesurfacemeshpoints were deEed bv local Bsulines. In the case of a uniform onedimensionh mesh, a 9spline with a displacement d centered at the mesh point i would produce displacements d / 4 at i + I and i  1 and zero elsewhere, while preserving COnhnuity of the first and second derivatives. Thus we can suppose that the discrete surfacedisplacement has the form
6S=Bd,
where B is a matrix with coefficients defined by the Bs lines, and di is the displacement associated with the
f s line centered at i. Then, using the discrete formulas, to f k t order the change in the cost is
61=BT6S=GTBd. Thus the gradient with respect to the Bspline coefficients is obtained by multiplying G by BT,and a descent step is defined by %tung
d=XBTG, 6S=Bd=XBBTG 61= J / G ( E , q ) aS(E,q) &dq where the gradient G 7) is Obtained bY evaluating the integralsin equation (57). Thus to reduce I we can choose 6S=XG
where A is suliiciently small and nonne ative In order to imuose a thickness constraint we can &fine a baseline surf&e So (6C) below which S (E, C) is not allowed to fall. Now we take X=X (4C) as a nonnegative function such that (58) S(S,C) + 6S(E,C) 2 So (t,0. Then the constraint is satisfied, while
where X is sufficiently smalland positive. The coefficients of B can be renormalizedto produce unit row sums. With a uniform mesh s acing in the computational domain this formula is uivafent to the use of a gradient modified by two passes3the explicit smoottung procedure
Withasimilarsmoothingpmcedureinthek discretization. Implicit smoothing may also be used. The smoothing equation
approximates the differential equation The costate solution (il is a legitimate test function for the weak form of the flow uations only if it is (tifferenuable. Smoothness *O be. PreFrvd the redesigned shape. It is therefore c~ciallymportant to mtmduce appropriate smoothmg procedures. In order to avoid disconunuitiesin the adjomt boundary condition which would be caused by the appearance of shock waves, the cost function for the target pressure may be mcdified
shea
Then
6I =
//
a ac ac
0  E+
If one sets 6S=Xc, then to first order the change in the cost
assuring an improvementif X IS sufficientl small andpositive, unless the process has already reacled a stationary
point at which G=O.
6.4 DesignofSwept WingsforVeryLowShockDrag
r
and the smootp quantity 2 replaces p  pd in the adjoint boundary conhhon.
The method has been used to cany out a stud of swept wing designs which might be a opriate for on range "osport aircraft. Since three Kensional calciations requm pbstantial computational resources, ,it is extremely important for $e practical implementatlon of the method to use fast soluhon algorithm for the flow and the
122
ndinini eouatinns. In this case the author's FL087 com__.__ .~ ~~~
~~~~
~
~~~
~
puter pro am has been used as the basis of the design method. KO87 solves the t&ee dimensional Euler equations with a cellcentered finite volume scheme, and uses residual averaeine and multierid acceleration to obtain ve rapid stedy k e solutio~s,~us,dly in25 to 50multigriycyc!es [66, 701.. Upwind biasing is used to produce nonoscillatory soluhons, and assure the clean capture of shock waves. This is introduced through the addition of carefully controlled numerical diffusion terms, with a magnitude of order Az3 in smooth parts of the flow. The adjoint equations are treated in the same way as the flow equations. The fluxes are first estimated by. central differences, and then modified b downwind biasing through numerical diffusive terms whch are supplied by the same subroutinesthat were used for the flow equatlons. The study has been focussed on wings designed for cruising at Mach 3 5 , with lift coefficients in the range of .5 to .55. In every case, the wing planform was fixed while !he sections were free to be chan ed arbitranly by the desi n method, with a restriction on lie minimum ttuckness. d e
rogram FL067. This program uses a cellvemx formufation, and has recently b e n mod!fied to. incorporate a local extremum dmnishmg a1 onthm with a ve level of numerical diffusion 174. m e ? run to fuVc'Z vergence it was found that a better estmate of the drag coefficient of the redesigned wing is 0.0094 at Mach 0.85 with a lift coefficient of 0.5, giving a lift to drag ratio of 53. The results from FL067 for the initial and final wings are illustrated in Figures 29 and 30. A calculation at Mach 0.500 shows a drag coefficient of 0.0087 for a lift coefficientof 0.5. Since in this case the flow is entirely subsonic, this provides an estimate of the vortex drag for this planform and lift distribution, which is just what one obtains from the standard formula for induced with an aspect ratio AR=9, and drag, CD=CL*/ETAR, an efficiency factor c=0.97. Thus the design method has reduced the shock wave drag coefficient to about 0.0007 at a lift coefficient of 0.5. Figure 31 shows the result of an analysis for an off design point with the Mach number increased to .86 with the same lift coefficient of .5. ? i s results in a flattopped pressure distribution terminahng with a weak shock of near1 uniform strength across the whole s an The drag coe&cient is ,0097. The penalty of ,000Lfis so small that this might be a preferred cruising condition. A second wing was designed in exactly ,the same manner as the first, starting from the same inmal geometry and with the same constraints, to give a l i i coefficientof .55 at Mach .85. This oroduces stroneer shock waves and is ~.~ therefore a more severe test ofthe hethod. In this case the total inviscid drag coefficient w+ reduced from0.0243 to 0.0134 in 40 design cycles. Agam the performance of the fmal design was verified by a calculation with FL067, and when theresult was fully converged the drag coefficient was found to be 0.0115. A subsonic calculation at Mach .500showsadragcoefficientofO.O107foraliiftcoefficient of 0.55. Thus in this case the shock wave drag coefficient is about 0.0008, For a representative transport aircrafl the parasite drag coefficient of the wing due to skin friction is about 0.0045. Also the fuselage drag coefficient is about 0.0050, the nacelle drag coefficient is about 0.0015. the empennage dra coefficient is about 0.0020.and excrescence drag coekcient is about 0.0010. This would give a total drag coefficient C~=0.0255for a l i i coefficient of 0.55. coresDondine to a lift to drag ratio LID=21.6. ' h s would be subsktial improvcmh over he values obtained by currently flying transport amraft. ~~~~~
0.6. This section,which has a thickness to chord ratio of 9.5 percent, was used at the ti Similar sections with an increased thickness were usefinboard. The variation of thickness was nonlinear with a more rapid increase near the rmt, where the thickness to chord raho of the basic section was multiplied by a factor of 1.47. The inboard sections were rotated upwards to give the initial wing 3.5 degrees twist fromroot to tip. ?e twodimensional pvssure distribution of the basic wing sectlon at its, desip int was introduced as a target pressufedistnbutlon uni& n l y across the span. This target is resumably not realizable, but serves to favor the estabhsK,ent of a relativel benign pressure distribution. The total inviscid dra coedcient, due to the combination of vortex and shoc wave drag, was also included in the cost function. Since the main objective of the study was to minimize the dra the @get pressure hstribution was reset after every fo design cycle to apispibuhon derived by smwthmg the exism pressure distnbuhon. Thm allows the scheme more freefom to make changes which reduce drag. The calculations were performed with the lift coefficient forced to approach a fixed value b adjusting the angle of attack every fifthiteration of the d w solution. It was found that the computational costs can be reduced b usin only 15 multigrid cycles in each flow solution, aniin ea& adjoint solution. Althou h this is not enough for full ,conv,ergence, it proves su%icient to provide a shape mdficahon wluch leads to an improvement. Figures 27 and 28 show a wing which was designed for a lift coefficient of .SO at Mach .85. In order to prevent the final wing from becoming too thin the threshold So (c, 7 ) ) was set at three quarters of the height of the bump S (t,7)) defining the initlal wing. This calculation was performed on a mesh with 192 intervals in the direction wrapping around the wine. 32 intervals in the normal n direction and 48 interval; in the spanwise C direction, giving a total of 2949 12 cells. The wing was specified b 33 scctions, each with 128 points, giving a total of 4234 design variables. The plots ihow ttie inifial wing geome pressure distribution, and the modified geometry an pres sure distributionafter 40 desi n cycles The total inviscid drag coefficient was reduced from 0.0210 to 0.01 12. The initlal design exhibits a very strong shock wave in the inboard region. It can be seen that this is corn letely eliminated, leaving a ve weak shock wave in &e outboard r r . TO verifyxe solution, the final geometry was an yzed with another method, using the computer
f
~~
r~~~~~~
~~~~
a
Optimization of Complex Configurations In order to treal more complex configurations one can use a numerical rid generation procedure to produce a bod fitted mesh or the initial geometry, and then modify mesh in sub uent design cycles by an analyhc perturbation formul3n the twodimensional case, for example, with computational coordinates c, q. let the boundary dish 6 m (0.Then the mesh olacement at n=O be 6 ~ (€1. points along the radial cobrdinateiinis €=constantcan be replaced by 6.5
F
de
yielding
7 "d
Such a procedure has been implemented by J. Reuther for the threedimensional Euler equations, and ap lied to the optimization of wingbody configurations[ 14!]. It is also possible to show that in the continuous limit the field integral in equation (57) can be eliminated. Let
123
the change in the coordinates 2, at k e d be 62i (5). Then, using the fact that the fluxes f, (w)satisfy the flow eauation (48). it is oossible to show bv a direct calculation
where A [email protected] derivation is given in,reference[781. Thus the perturbahon equahon can be wntten as
be significantly improved by innovative concepts, such as the idea of time inclining. It can be anticipated that interdisciplinary applications in which CFD ii cou led with the com utauonal analysis of other properties ofthe desien will D av an increasinelv imoortant role. These appkations'miy include sm%ral, thermal and electromagnetic analysis. Aeroelastic problems and integrated control system and aerodynamic design are likely target areas. The development of improvedalgorithmsconh
P
~~~~~~~~
~~~~
~~~~~
.
I
where 6w is the variation in the solution at fixed caused by the change in the boundary, whle 6w' is the change in the original solution w(t) corresponding to the mesh movement 62 (t)

I 
Figure 19: Concept for a numerical wind tunnel. Now
and if w satisfies the adjoint e uation the entire field integral is eliminated, leaving on?y the boundary integral in equation (57). In an actual discretization the field terms are not zero. but this result suggests that they should be small if a fine enough mesh is used, and mi ht be dropped. This allows a drastic simplification ofthe treatment of complex con6 urations. Pieliminary numerical experiments with airfoi and wing calculations indicate roughly the same convereence wlth and without the field terms in the eradient.
Figure 2 0 Advanced numerical wind tunnel.
ues to be impr!ant in providing the basic building blocks for numerical simulauon. In narticular. better ...... error ..... esti..mation procedures must be develo d and incorporated in the smulation software to provi error control. The basic simulation software is only one of the needed ingredients, however. The flow solviz must be embedded in a userfriendly system for geometry modelin output analysis, and data management that will provife a corn lete 7. OUTLOOK AND CONCLUSIONS ients numerical desi n environment. These are the ingre$ which are n d e d for the fuy realization of the concept of Betteralgorithmsandbettercomputerhardwarehavecon a numerical wind tunnel. Figures 19 and 20 illustrate the tributed about equally to the progress of computational way in which a numencal wind tunnel might evolve from science in the last two decades. In 1970 the Control Data current techniques, which involve massive data handling 6600 represented the state of the art in computer hardtasks,to a fully integrated design envlronment. ware with a s ed of about IO6 operations per second In the long run, corn utational simulation should become (one megafloprwhile in 1990 the 8 processor Cray YMP the pMci~altoo! o ,the aerodynamic design,processbeoffered a performance of about lo9 operations per seccause of e flehibiI!ty it provides for the rapid and comond (one gigaflop). Correspondingly, steadystate Euler paratively inexpensive evaluation of alternative designs. calculations which resuired 5.OOOlO.OOO stem nrior to ind because it can be inrefated in anumerical designen1980 could be perfo&ed in IC50 steps in i99b using vironmentpviding for 0th mulu dsciphnary analysis mulugrid acceleration. With the advent of massively parand multi isciplinary optimization. allel computers it appears that the progress of computer hardware may even accelerate. Teraflop machines offering further iniprovement by a factor of I ,OOO are likely to REFERENCES be available within a few years. Parallel architectures will force a reappraisal of existin algorithms, and their effecIll R. Abid, C.G. Speziale, and S. Thangam. Aptive utilizauon will require extensive development of lication of a new kr model to near wall turnew parallel software. gulent flows. AIAA Paper 910614, AIAA 29th Aerospace Sciences Meeting, Reno, NV, January In arallel with the transition to more sophisticated algo1991. ri&s, the present challengeis to extend the effective use of CFJJ to more complex applications. A key roblem is 121 H. Aiso. Admissibility of difference a roxima the treatment of multi le s ace and time scapes. These tions for scalar conservauon laws. f % r o s h i ~ anse not onlv in turbu?ent Eows. but also in manv other Math. Journal, 23:1561,1993. situations suih as chemically reactin flows, com6ustion. flame fronts and lasma dynamics. k o t h e r challenge, is [3] J. J. Alonso and A.,Jameson. Full implicit timepresented by proilems with moving boundanes. Eh?marchin aeroelashc so~utions. pa er 94ples include helicopter rotors. and rotorstator interacuon 0056, Af4A 32nd Aems ace Sciences hfeeting, in turbomachinery: Algorithms for these problems can Reno, Nevada, January 1954. ~~~~~~~
~~~~~~~~~~
~~~
%

I
P
8e
&M
124
[4] J. J. Alonso, L. Martinelli, and A. Jameson. Mu!tigrid unsteady NavierStokes calculauons with aeroelastic a phcations. AlAA aper 95 0048, AIAA ~ 3 r dlerospace Sciences Geeting, Reno, Nevada, January 1995.
[ZO] J.A. Benek, T.L. Donegan, apd N.E. Sub+ Extended Chunera gnd embeddm scheme with applications to viscous flows. AI.& Paper 871 126, AIAA 8th Computational Fluid Dynamics Conference. Honolulu, HI, 1987.
[5] B.K. Anderson, J.L.Thomas,,+d B. Van Leer. A companson of flux vector sphtUngs for the Euler equations. AIAA Paper 850122, Reno, NV,January 1985.
[211 M. Ber er and A. Jameson. Automatic ada tive grid reAement for the Euler equations.AIAA !ournul, 23561568,1985.
[6] W.K.Anderson, J.L. Thomas, and D,L. Whiffield. Multigrid accelerationof the flux s lit Euler equations.AIAA Paper 860274,AIAA 4th Aerospace Sciences Meeting, Reno, January 1986.
s
[7] T.J. Baker. Mesh generation by a sequence of transformations. Appl. Num. Math., 2515528,1986. [8] N. Balakrishnan and S . M. Deshpande. New upwind schemes wjth waveparUcle s littin for in12, Inviscid compressible flows. Reporr 81 dian Institute of Science, 1991.
Id
a
[9] B. Baldwin and H. L o a . Thin layer a proximation and algebraic model for separate turbulent flow. AIAA Paper 78257.1978. [IO] B.S. Baldwin and T.J.Barth. A onee uation turbulence trans rt model for high ReynJds number wallboundeg0flows. AIAA Paper 910610, AIAA 29th Aerospace SciencesMeeung, Reno,NV,January 1991. [ l l ] T. J. Barth. Aspects of unstructured ds and finite volume solversfor the Euler and avier Stokes equations. In von Karman Institutefor namics Lecture Series Notes 199405, 1994.
P.
[ 121 T.J. Barth and P.O. Frederickson. Higher order so
lution of the Euler equauons on unstructured gnds uadratic reconstrucuon. AIAA paper 90anuary 1990.
E3 3
The design [13] T.J. Barth and D.C. Jes rsen and application of upwmrxhemes on unstructured meshes. AIAA paper 890366, AIAA 27th Aeros ace Sciences Meeting, Reno,Nevada, January 1689. [14] J.T. Batina. Implicit, fluxsplit Euler schemes for unsteady aerodynarmc analysis mvolvln unstructureddynarmcmeshes. AIAApaper900836.ApnI 1990. [ 151 E Bauer, P. Garakdian, D. Kom, and A. Jameson.
Supercritical Wng Sections 11. Spnnger Verlag, New York, 1975.
[16] 0. Baysal. e d M. E.,Eleshaky. , Aerodynamic design optmuanon usmg sensiuvi anaysis and corn utational fluid dynamics. ' paper 91047f 29th Aeros ace Sciences Meeung, Reno, Nevada, January 1891. 1171 ., R.W. Beam and R.F. Warmine. An imdicit finite differencealgorithmforhyer6oIic s scams in conservauon form. J. Comp. Bhys., 23:871 IO, 1976.
'A
..
~~~
~~~
Reno, Nevada, January 1995. [19] J.A. Benek, P.G. Bunin and J.L. Steger. A 3D Chimera d embed?& techni ue In Pmceedings 7th ComputarioMI h ~ i Dynamh ics Conference, pages 507512, Cincinnati. OH, 1985. AIAAPaper 851523.
d
1221. M. BergerandRJ. Levewe. Anadaptivecartesian . mesh algorithm for the Euler e uations in arbitrary geometnes. AIAA Paper 891830.1989.
[23] J.P. Boris and D.L. Book. Flux corrected transport, 1 SHASTA, a fluid trans ort al orithm that works. J. Comp. fhys., I1:38& 1975. [24] A. Brandt. Multileveladaptivesolutionsto boundvalue problems. Math. Comp., 31:333390, %7. M.O. Bristeau, R. Glowinski, J. Periaux. P. Pemer, 0. Pironneau, and C. Pokier. On the numencal solution of nonlinear roblems in fluid d y n m c s by least squares and Knit, element methods (U), application applicatioi to transonic flow simulations. Comp. comp. Meth Appl. Mech. andEng., 51:363394,1985. [26] R.J. Busch, Jr. Corn utational fluid dynamics in the desi n of the No&o /McDonnell Dou las YF23 prototype., paper 911629, AIAA 21st Flud D narmcs, Plasmgynamics & Lasers Conference. Aonolulu. Hawau. 1991.
h ,F
A h
8.
[271 C. Canuto, M.Y. Hussaini, A. uarteroni, and D.A. Zan Spectral Methods in urd DyMmics. Springer+erIag, 1987. [28] M.H. Carpenter, D. Gottlieb, ,and S . Abarbanel. Timestable boundary conditions for f i ~ t e difference schemes solving hyperbolic systems: Methodology and a lication to hi h order comact schemes. C tA Ir& oR pe! 939, hmpton, VA, barch 1993. 1291 D.A. Caughey. A diagonal implicit multigrid algorithmfor the Euler equations. AIAA Paper 87453, 25th Aerospace Sciences Meeting, Reno, January 1987.
f
[30] T. Cebeci and A.M.O. Smith. Analysis o Turbulent Boundary Layers. Acadenuc Press, 19 4. 1311 S.R. Chmvarthy.,Relaxation methods for unfactored impkit upwind schemes. AIAA Paper 840165, AIAA 22nd Aerospace Sciences Meeting, Reno,January 1984. [32] S.R. Chakravarth A Harten, and S., Osher. Essentially nonoscaatory shock of uniformly very hi h 0339, AIAA 24th bferospace Reno, January 1986. [33] R. Chipman and A. Jameson. Fully conservative numencal solutions for unsteady irrotational transonic flow about airfoils. AIAA Paper 791555, AIAA 12th Fluid and Plasma D namics Conference, W a r n b u r g , VA, July 19&. 1341 S.E. C W and S.D.,Thomas. Euler/experimental correlations of SONC boom ressure signatures. AIAA Paper 913276. AIAA 8 th Applied Aerodynamics Conference, Baltimore, September 1991. [351 T.J. Coakley. Numerical simulation of viscous transonic airfoil flows. AIAA Paper 870416. AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.
125
[36] J.P. CroisiLle and P. Villedieu. Kinetic flux spli,tting schemes for hypersomc flows. In M. Napobtano and F. Sobetta, editors, Pmc 13th International Congress on Numerical Methods in FluidDynamics, Dazes 31C3313. Rome, July 1992. Springer [37] R.M. Cummings,Y.M.Rirk,L.B. Schiff,andN.M. Chadeqian. NavierStokes,predictionsfor the F 18 wing and fuselage a1 large  incidence. J. ofAircraj?, 29:<65574.1992. [38] G.Dahlquist. A s cial stability roblemforlinear multistep m e t h z B l T , 3:274! 1963. [39] J.F. Dannenhoffer and J.R. Baron. Robust adautation for comvlex transonic flows. Papk 860495, A L h 24th Aerospace Sciences Meeting, Reno,January 1986.
[53] W.Hackbusch. On the multigrid method a plied to hfference equations. Compuring, 2029f306, 1978. [54] M.Hafez, J:C. South,andE.M. Man. Artificial compressibility method for numencal soluuons of the transonicfull otentialequation.AIAA Journal. 17:838844,197!. [55] M.G. Hall. Cell vertex multi d schemes for solution of the Euler uauons. E P m c . IMA Conference on Numeric3 Methods for Fluid Dynamics, Reading, Apnl 1985.
ad
[56] A. Harten. High resolution schemes for hyperbolic conservation laws. J. Comp. Phys.. 49357393, 1983.
[40]D. Deganimd and L. Schiff. Computation Of ~ U I ~ U 
[57] P.W. Hemker and S.P.Spekregse. Mulugrid solution of the steady Euler equatlons. In Pmc. Oberwol ach Meeting on Multigrid Methods. December 19d.
[41] B. Delaunay. Sur la sphhre vide. Bull. Acad. Science USSR VII: Class Scil, Mat. Nat., pages 793800.1934.
[58] J.L. Hess and A.M.O. Smith. Cgculation of nonWing otential flow abou! arbitrary three dimensionafbcdies. Douglas Arcraft Report ES 40622,1962.
lent supersom flows around pointed bodies havin crossflowseparation. 1. Comp. Phys.. 6617319t 1986.
[42] S.M. Deshpande. On the Maxwellian distribution, symmetnc form and entropy conservation for the Euler equations. NASA TP 2583,1986. [43] A. Eberle. A finite volume method for calculating transonic potential flow around wings from the minimum pressure inte al. NASA TM 75324, 1978. Translated from &B UFE 1407(0). [44] P.R. Eiseman., A multisurface method of coordinate generauon. 1. Comp. Phys., 33:118150, 1979.
[46] L.E. Eriksson. Generatjon of boundaryconforming gnds ,around wingbod configurations usin transfimte mterpolauon. AA Journal, 201313k320,1982.
Al
namics Conference, San Diego, CA, June 1995. [48] R.P. Fedorenko. The s d of conver ence of one Matf and Math. iterative process. US&?Comp. Phys., 4:227235,1964. [49] C.W. Gear. The numerical integration of stiff ordin a y differentialequations. Report 221,. Universi of h o i s Department of Computer Science, 1967 [501 M. Giles, M. Drela. and W.T. Thompkins. Newton solution of direct and inverse transonic Euler equations. AIAA Paper 851530, Cincinnati, 1985.
[51] S.K.Godunov. A difference method for the numerical calculation of discontinuous solutions of hydrodynamic equations. Mat. Sbomik, 47271306,1959. Translated as P R S 7225 by U.S. Dept. of Commerce, 1960. [52] M.H. Ha. The im act of turbulence modelling on the numerical prJctionof flows. In M. Na olitano and F. Solbetta, editors, Pmc. o the 13th fnternational Conference on Numerica Methods in Fluid DyFmics, pages 2746, Rome, Italy, July 1992. Spnnger Verlag, 1993.
f
[59] C. Hirsch, C.Lac01, andH. Deconinck. Convection algorithms based on a dia onahation procedure forthernultidimensional~~re uauons ~n ~mc. AIM 8th CO utational Flu!d hnamics C a ence. Dazes 3  6 7 6 , Hawau, June 1987. Paper ST1163. [60] C. Hirsch and P. Van Ransbmk. Multidimensional upwindin and artificial dissi ation. Technical re ort Publfshed in Fmntiers o Computational duid'Dynarnics 1994, D.A. J u g h e y and M. M. Hafez, editors, Wiley, pp. 597626. [61] D.G. Holmes and S.H.Lapon. Ada tive Irianular meshes for compressible flow sofutions. In hvceedings First International Conferenceon Numerical Grid Generation in Corn utational Fluid namics, pages 413424, Landhut, FRG,July 386. [62] T.J.R. Hughes, L.P. Frapca, andM. Mallet. A new finite element formulauon for computauonal fluid d namics.1, Symmetric forms of the compressible d l e r and NavierStokes equations and the second law of thermodynamics. Comp. Meth. Appl. Mech and Eng., 59:223231,1986. [63] A. Jameson. Iterative solution of transonic flows over airfoils and wings, including flows at Mach 1. C o r n on Pure and Appl. Math., 27:283309, 1974.
[64]A. Jameson. Transonic potential flow calculauons in conservation form. In Pmc. AIAA 2nd Com U rational ~ ~ uDynamics i d Conference, pages I& 161. Hartford, 1975. [65] A. Jameson. Solution of the Euler equations by a mulugnd method. Appl. Math. and Comp., 13:327356.1983. [66] A. Jameson. Solutionof the Euler e uations for two dimensional transon~cflow by a mitigrid method. Appl. Math Comp., 13:327356,1983. [67l A. Jameson. Multigrid al orithms for compressible flow calculations. In % e c o d Eumpean Conference on Multigrid Methods, Cologne, October 1985. Princeton University Report MAE 1743.
126
I681 A. Jameson. Nonoscillatory shock c turin scheme using flux limited dissi ation. B.2 Engquist, S . Osher, and R.C.J. %mmerville, editors, Lectures in Applied Mathematics, Vol. 22, Pari I , Lorge Scale Computations in Fluid Mechanics, pages 345370.AMs, 1985.
%
[821 A. Jameson,,T.J.Baker,apdN.P.Weatherill. Calculation of inviscid transonic flow over acomplete aircraft. AIAAPaper860103,AIAA24thAerospace Sciences Meeung, Reno, January 1986. [83] A. Jameson and D.J. Mavriplis. Multigrid solution of the Euler eauations on unstructured and adantive 'ds. In S . MZCormick, editor, MultigridMethods, Kay,A lications and~upercomputing.~ecture Notes in Ere and Ap lied Mathematics, volume I IO,pages 413430,Af)pril1987. [841 A. Jameson, ,W. Schmidt, and E. 'hrkel. Numencal solutlon of the Euler equations by finite volume methods using RungeKutta time stepping schemes. AIAA Paper811259,1981. ~
[69] A. Jameson. Transonic flow calculations for aircraft. In F. Brezzi. editor, Lecture Notes in Marhematics, Numerical Methods in Fluid Dynamics, pages 156242. SpringerVerlag, 1985. [70]A. Jameson. Multi 'd al onthmsforcompressible flow calculations. F W . 8ackbuschand U. Trottenber , editors, Lecture Notes in Mathematics, Vol. I 2 h ~ g e 166201. s PrFeedin softhe2ndEuropean onference on Mulugrid detbods, Cologne, 1985,SpringerVerlag, 1986. [71]A. Jameson. A vertex based multigrid algorithm for threedimensional corn ressible flow calculations. In T.E. Tezduar and F.J.R. Hu hes, editors Numerical Methodr for Cam ressib& Flow  Fi: nite Di erence, Element Al k tVolume Techniques, 1986. S M E Publication AMD 78.
f
1721 A. Jameson. Aercd namic desi n via control theory. J. Sci. Comp., $233260, b88. [73] A. Jameson. Automatic design of mansonic airfoils to reduce the shock induced pressure drag. In Proceedings of the 3Ist Ismel Annual Conference on Aviation and Aemnautics, Tel Aviv, pages 517, February 1990. [74] A. Jameson. lime dependent calculations using multi 'd, with applicatlons to unsteady flows ast airfoipand wings. AIAA aper 91 1596, loth Corn utational Fluid gynamics Conference, Honolulu.!hawaii. June 1991.
A%A
[75]A. Jameson. Artificial diffusion, upwind biFing, limiters and their effect on accuracy and mulugnd convergence in transonic and hypersonic flows. AIAA Paper 933359,AIAA 11th Computational Fluid Dynamics Conference, Orlando, FL. July 1993. [76] A. Jameson. Artificial diffusion. upwind biasing, limiters and their effect on accuracy and,mulud convergence in transomc and hy BONC flow. paper 933359. AIAA I ~ t gmputatiqnal h Fluid D namics Conference, Orlando, Flonda, July 1993.
.EM
[77]A. Jameson. Optimum aerod namic desi n via boundary control. Technicdre ort, A6ARp FDPNon K m a n hstiFtespecial 8,urseon opumum Desi n Methods m Aerodynmcs, Brussels, April 199j [781 A. Jameson. MAE Technical Report 2050, F'rinceton University, F'rinceton, New Jersey, October 1995. [791 A. Jameson. Analysis and design of numerical schemes for, gas d y n m c s 1, e f i c i a l diffusion, upwind biasmg, b t e r s and theu effect on multigrid convergence. Int. J. of Comp. Fluid Dyn., 4171218,1995. [80]A. Jameson. Analysis and design of numerical schemes for as dynamics 2, artificial diffusion and dscrete dock shuchue. Int. J. of Comp. Fluid Dyn.. To Appear. 1811 A. Jameson and T.J. Baker. Improvements to the aircraft Euler method. AIAA Paper 870452, AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.
~~
~
~
~
~~~
~~~~
~~
~
1851 . . A. Jameson. W. Schmidt. andE. Turkel. Numerical
solutions of the Euler equations by finite volume methods with Run e Kutta time ste ping schemes. AIAA paper 8I12%, January I98!I
[86] A. Jameson and E. Turkel. Implicit schemes and LU decompositions. Math. Comp., 37385397. 1981. [87]M. Ja a r m and A. Jameson. Multigrid solution of the dvierStokes uations for flow over wings. AIAA aper 88 073,AlAA 26th Aeros ace Sci encesbeeting,Reno, Nevada, January 1888. [881 D. Johnson and L. King. A mathematically simple turbulence closure model for attached and separated turbulent boundary layers. A A 4 Journal, 23:16841692,1985. [891 W.P. Jones and B.E. Launder. The calculation of lowReynoldsnumber phenomena with a twoe uationmodelof turbulence. Int. J. ofHeat Tron., 1%:1119l130,1973. [90]W.H. Jou. Bwin Memorandum AEROB113BL92018,Septemfer 1992. To Joseph Shang. [91] T.J. Kao, T.Y. Su,andN.J. Yu. Navierstokescalculations for transport wingbody configurations with nacelles and shuts. AIAA Paper 932945, AIAA 24th Fluid Dynamics Conference, Orlando, July 1993. [921 G.E. Karniadakis and S.A. Orszag. Nodes, modes and flow codes. Physics Today, pages 3442, March 1993. [931 M.H. Lallemand and A. Dervieux. A multi'd finiteelement method for solvin the two%ensional Euler equations. In S.F. &Connick, editor, Pmceedin s of the Third Cop er Mountain Con erence on &It, rid Methods, &care Notes in ure and AppJiedhathematics, pages 337363, Copper Mountam. Apnl1987. [94]A.M. Landsber ,J P Boris. W. Sandberg, and T.R. Young. Naval sh$ ;uperstructure desi threedimension flows using an le1 method. High Perfomnce Com uting ?99$ Gmnd Challenges in ComputerSimu%tion, 1993.
d
e&i2tmgF
[95]P. D. Lax. H perbolic systems of conservation laws. SIAM dgional Series on Appl. Math., U, 1973.
[97]P.D. LaxandB. Wendroff. Systems ofconservation laws. Comm Pure. Appl. Math., 13917237,1960. [98] B. Van Leer. Towards the ultimate conservative differencescheme. U. Monotonicit and conservation combined in a second order scgeme. J. Comp. Phys., 14:361370,1974,
127
[99] B. Van Leer. Towards the ultimate conservative difference scheme. ID upstreamcentered finitedifference schemes for ideal compressible flow. J. Comp. Phys., 23963275,1975.
[ 1011 B. Van Leer. P r o ~ s &multi+imension_alups wind differencing. In M. Napolltano and F. Sol
betta, editors, Pmc. 13th International Conference on Numerical Methotds in Fluid D f u z m y jages 126, Rome, July 1992. !jpnnger erlag 19 3
[lo21 B. Van Leer, W. T. Lee, and P. L. Roe. Characteristic tlme stepping or local preconditionin of the Euler equauons. AIAA aper 91 1552, A h A loth Corn utational Fluid &namicE Conference, Honolulu,%awaii, June 1991. [ 1031 D. Lefebre, J.Peraire, and K.Mor an. Fhte ele
ment least squares solutions of the, E! uler equations using linear and quadrauc ap roxunauons. Int. J. Comp. FluidDynamics, 1:l$3. 1993.
[ 1151 L. Mdnelli, A. Jameson, +d E. Malfa. Numerical
simulation of threedimennonal vortex flows over delta wing configurations. In M. Napolitano and E Solbetta, editors, Pmc. 13th International Confrence on Numerical Methoh in Fluid Dynamics, ages 534538. Rome, Italy, July 1992. Springer eerlag, 1993.
[116] L. Martinelli and V. Yakhot. Y G  b F e d ,turbulence transport approximations w~thapphcauons to transonic flows. AIAA Paper 891950, AIAA 9th Com utational Fluid Dynamics Conference, Buff a l o , b , June 1989. [I171 D.J. Mavri lis and A. Jameson, Multifid solution of the%avierStokes equauons on tnangular meshes. AIAA Journal, 28(8):14151425, August 1990. [118] D.J. Mawiplis and L. Martinelli. Multigrid solution of compressible turbulentflow on unstructured meshes using a two uation model. AIAA Paper 910237, January 1 9 3 . [119] N. D. Melson, M. D. Sanehilt, and H. L. Atkins. limeaccurate NavierStokes calculauons w~th multigrid acceleration. In Pmceedings of the Sixth Copper Mountain Conference on Multigrid Methods, Copper Mountam, Apnl1993.
[ 1041 S.K.Lele. Compact finite difference schemes with
trallike resoluuon. ?KO.
CTR Manuscnpt 107,
[lo51 J.L. Lions.
[106] MS. Liou and C.J. Steffen. A new flux splitting scheme. J. Comp. Phys., 107:2339,1993.
[lo7 E Liu and A. Jameson. Mu!tigrid NavierStokes calculationsfor threedimensional cascades. AlAA aper 920190, AIAA 30th Aeros ace Sciences beeting, Reno, Nevada, January 1982.
[121] E Menter. Zonal twc uation k w turbulence models for aerod namicYows. ~ I A A Pcper 932906, AIAA 2 4 d Fluid D y n m c s Meehng, Orlando, July 1993. [ 1221 E.M. Murman. Analysis of embedded shock waves
calculated by relaxation methods. AIM Journal, 12626433.1974.
[ 1231 E.M. Murman and J.D. Cole. Calculation of plane
steady transonic flows. AIAA Journal, 9:114121. 1971.
[ 1081 R. Lohner and D. Marfin. An implicit lineletbased
solver for incompressible flows. AIAA a er 92 0668, AIAA 30th Aerospace Sciences?&etinL Reno,W.January 1992.
[ 1241 R.H. Ni. A multiple grid scheme for solving the Eu
lerequations. AlAA Journal,2015651571,1982.
[lo91 R. Lohner, K. Morgan, and,]. Peraire. Improved adaptive rehement strateges for the finite element aerod namic conligurations. AIAA Paper 860499, Al%A 24th Aerospace Sciences Meeting, Reno, January 1986.
[ 1251
[I101 R. Lohner, K. Morgan, J. Peraire, and O.C. Zienkiewicz. F~mteelement methods for high s ed flows. In Pmc. ALAA 7th Computational Xid Dynamics conference. Cincinnati, OH, 1985. AIAAPaper851531.
[126] J.T. Oden. L. Demkowicz, T. Liszka, and W.Rachowicz. h ad tive finite e!ement methods for compressibye an? incompressibleflows. In S. L. Venneri A. K.Noor, editor, Pmceedings of the Symposium on Computational Technology on Flight Vehicles, a es 523534, Washington, D.C., November l d & r g a m o n .
[l111 R. Liihner and P. Parikh. Generation of three.dimensional unstructured grids by the advancing front method. AIAA Paper 880515. Reno, W , January 1988.
S. 0bayas.F and K. Kuwakara. LU factorization
of an imphcit scheme for the compressible NavierStokes uations. AIAA Paper 841670, AIAA 17th F l 8 Dynamics and Plasma Dynamics Conference, Snowmass.June 1984.
[127] S. Orszag and D. Gottlieh. Numerical analysis of spectral methods. SIAM Regional Series on Appl. Math., 26,1977.
[ 1121 R.W. MacCormack ,and A.J. ,Paullay. . Com uta
11281 S. Osher. , Riemann solvers, fhe entro
[ I 131 A. Majda and S: Osher. Numerical viscosity and theentro condrtlon. C o r n onPureAppl. Math., 32:797& 1979.
[129] S. Osher and S. Chakravarthy., Hi h resolution schemes and the entropy condmon. h U J . Num Anal., 21:955984,1984.
[114] L.Martine1liandA. Jameson.Validationofamultid method for the Reynolds averaged equauons. E 4 A paper880414.1988.
[130] S. Osher and E Solomon. Upwind difference schemes for hyperbolic s stems of conservation laws. Math. Comp., 38:336374, 1982.
tionalefiiciency acheved by m e sphmng of E n i i differenceoperators.AIAAPaper 72154.1972.
condition, and difference approxunatlons.S l A r J . Num. Anal., 121:217235,1984.
128
U1311 S. Osher and E. Tadmor. On the convergence of difference approximations to scalar conservation laws. Math. Comp., 501951, 1988. [I321 H. Pailhe and H. Deconinck. A review of multidimensional upwind residual distribution schemes for the euler equations. To appear in CFD Review, 1995. [1331 N. A. Pierce and M. B. Giles. Preconditioning on stretchedmeshes. Report9YlO 1995, OxfordUniversity Computing Laboratory, Oxford, December 1995. U341 0. Pironneau. Optimal Sha e Desi n or ENiptic System. SpnngerVerlag, Yorf, f984.
dw
[ 1351 K.G. Powell and B. van Leer. A genuinely multidi
[1481 C.L. Rumsey andYN. Vatsa. A corn arison of the predictive capabilities of several tur ulence models using upwind and centered  difference computer codes. AIAA Paper 930192, AIAA 31st Aerospace Sciences Meeting, Reno, January 1993.
%
11491 S.S. Samant, J.E. Bussoletti, ET. Johnson, R.H. Burkhart,B.L. Everson, R.G. Melvin, D.P. Youn L.L. Erickson, and M.D. Madson. TRANAIR corn uter code for transonic anal ses of a r b i t r ~ contf)gurations.AIAA Paper 870834.1987,
k
[ 1501 K. Sawada and S. Takanashi. A numerical investi
ation on wing/nacelle interferences of USB conEgFation. In Pmceedin s AIAA 25th Aemspace Sciences Meeting, Reno,&, 1987. AIAA paper 870455.
mensional upwind cellvertex scheme for the Euler equations. AIAA Paper 890095, AIAA 27th Aerospace SciencesMeeting, Reno, January 1989.
[151] C.W. Shu and S. Osher. Efficient implementation of essentially nonoscillato shockca turing schemes. J. Comp. Phys., 77:43v471,198!.
[136] K. Prendergast and K. Xu. Numerical hydrodynamics from gas kmetic theory. J. Comp. Phys., 109:5366, November 1993. [ 1371 T.H. pulliam and J.L. Steger. Implicit finite difference simulations of threedimensional compressible flow. AlAA Journal, 18:159167,1980,
[1521 C.W. Shu and S. Osher. Efficient implementation of essentially nonoscillato shockca turing schemes 11. J. Comp. Phys., 83:%78,1988.
[I381 J.J. Quirk. An altemauve to unstructured gnds for compuung as dynamics flowsnbout arbitranly com lex twocknensional bodies. ICASE Report 92$Hampton. VA, February 1992. 11391 R. Radespiel, C. Rossow, and R.C. Swanson. An efficient cellvertex multigrid scheme for the threedimensional NavierStokes equations. In Pmc. AIAA 9th Corn utational Fluid D namics Conference, pages %9260, Buffalo, I$Y, June 1989. AIAA Paper 891953CP. U1401 M.M. RaiandP. Moin. Direct numerical simulation of transition and turbulence in a spatially evolvin boundary layer. AIAA Pa er 91 1607 CP, loth Corn utational Fluidgynakcs Conference, Honolulu,&. June 1991.
d
[1531 C.W. Sbu, T.A. Zang. G. Er1ebacher.D. Whitaker, and S. Osher. Highorder EN0 schemes ap lied to two and threedimensional compressible iow, Appl. Num Math., 9:4571, 1992. [I541 B. R Smith. A near wall model for the IC  1 two equation turbulence model. AIAA paper 942386, 25th AIAA Fluid D namics Conference, Colorado Springs, CO, June 1'994. [I551 L.M. Smith and W.C. Reynolds. On the YakbotOrszag renormalization ou for deriving turbu lence statistics and modef PRys. FluidSA, 4:36L 390,1992. [1561 R.E. Smith. Threedimensional aleebraic mesh eneration. In Pmc. AlAA 6th Corn uktional Fluid bynamics Con erence, Danvers, d A , 1983. AIAA Paper 8 3  1 9 d [ 1571 R.L. Sorenson. Elliptic generauon of compressible
threedimensional gnds about realistic aircraft. In J. HauserandC.Taylor, editors,Infemotiono/Con
(1411 S.V.RaoandS.M.Deshpande. Aclassofefficient kinetic U wind methods for com ressible flows. Report9fFMI1, IndianInstituteo! Science, 1991.
feynce on Numerical Grid Generation in Com U fafionalF/uidDynomics,Landshut, F. R. G., 1886
[ 1421 J. Reuther and A. Jameson. Aerod namic shape op
[ 1581 R.L. Sorenson. Tbreedimensional elliptic gjid
timization of wing and wingbod; c o d using control theory. AIAA paper 9 5  0 E % : Aerospace Sciences Meeting and Exibit, Reno, Nevada, January 1995.
H431 J. Reutherand A. Jameson. Aerod namic shape optimization of wing and wingbo y c o d urations usin control theory. AIAA aper 9501f7, AIAA 331Lf Aerospace Sciences d e t i n g , Reno,Nevada, January 1995.
J
[1441 H. Rieger q d A. Jameson: Solution of steady threehmensional comgressible Euler and NavierStokes e uations b an im licit LU scheme. AIAA p e r ' 8 h 6 1 9 , A h 26% Aeros ace Sciences ceeting, Reno, Nevada, January 19118. 11451 A. Rizzi and L.E. Eriksson. Computation of flow around wings based on theEulerequations.J. Fluid Mech., 148:4571,1984. [146] P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comp. Phys., 43:357372.198 1. [I471 P.E. Rubbert and G.R. Saaris. A general threedimensional otential flow method a lied to V/STOL aeroiynamics. SAE Paper 680&, 1968.
eneration for an F16. In J.L. Steger and B.Generation F. Thompson, editors, ThneDimensiona/ ~ r i d or Corn lex Confi umrions: Recent Progress, lh8.AGhiDograpf?
[159] P. Spalart and S. Allmaras. A oneequation turbulent model for aerodynamic flows. AIAA Paper 920439, AIAA 30th Aerospace SciencesMeeting, Reno, Nv, January 1992. [160] C.G. Speziale,, E.C. Anderson, and R. Abid. A critical evaluatlon of twoequatlon models for near wall turbulence. AIAA Paper9C1481, June 1990. ICASE Report 9046. [I611 J.L. StegerandDS Chaussee. Generation ofbodv.fitted coordinates usin h perbolic artial diffirential e uations. S I A l ~s c i . ank'srat. comp., 1:4314\7, 1980. ~~
~
[162] J.L. StegerandR.F. Warmine. Flux vectorsnlitting of the iniisid as dynamic equations with ap =:= hca tions to finite Lfifferencemethods. J. Comp. Jhys.: 40:263293.1981. Cl631 B. Stoufflet,J. Periaux, F. Fezoui, and A. Dervieux. Numerical simulation of 3D hypersonic Euler flows around space vehicles using adapted finite elements. AIAA paper 870560, January 1987.
129
I1641 G.StrangandG.Fix.AnaI sisoftheFiniteEIement Method. Prentice Hall, 1 6 3 . 11653 R.C. Swanson and E. Turkel. On centraldiffercnce and U wind schemes. J. Comp. Phys., 101:297306. P992. [166] P.K. Sweby. High resolution schemes usin flux limiters for hyperbolicconservation laws. S h J. NumAnaL,21:9951011.1984. U671 S. Ta’asan, G.Kuruvila. and M. D. Salas. Aerod namic desi n and optimization in one shot. paper 92&, 30th Aeros ace Sciences Meeting and Enbit, Reno, Nevada, %nuary 1992.
AIL
11681 Shlomo Ta’asan. canonical foms of multidimensional, steady inviscid flows. ICASE 9334. Institute for Computer A lications in Science and Engineering, Hampton,vA, 1993. [169] E. T+,mor. Numerical viscosi and the entropy condmon for conservauve erence schemes. Math. Comp., 32369382,1984.
d2
[170] J.F. Ihompson, F.C. Thames,,and C.W. Mastin. Automatic numerical generation of Myfitted curvilinear coordinate system for field containing any number of arbitrary twodimensional bodies. J. Comp. Phys., 15999319,1974. 11711 J.F. .~
Ihomnson. Z.U.A. Warsi. and C.W. Mastin. Boundaryhted coordinate s stems for numerical solution of artial differentide uations: A review. J. Comp. fiys.. 47:l108. 1989.
[ 1721 E. Turkel. Preconditioned methods for solving the
incom ressible and low speed equations. J. Comp. Phys.,!J2:277298.1987.
[173] V. Venkatakrishnan. Newton solution of inviscid and viscous problems. AIAA paper 880413, January 1988. I1741 V. Venkatakrishnan. Convergence to steady state solutions of the Euler uatlons on unstructured ds with limiters, A l a aper 930880, AIAA %st Aerospace Sciences deting,  Reno, Nevada, January 1993. [175] G.Voronoi. Nouvelles plications des parametres continus a la theoriexs formes quadrati ues Deuxihe &moire: Recherches sur les ar%el: loedres rimitifs. J. Reine Angav. Math., 134:198287,
, a8.
[176] Y.Wada and MS. Liou. A flux splitting scheme withhi h resolution androbustness fordmontinuities. &$aperWOO83. AIAA 32ndAeros ace Sciences eehng. Reno,Nevada, January 1984. [177] N.P. Weatherill and C.A. Forsey. Grid generation and flow calculations for aircraft geometries. J. Aircraft, 22855860.1985. [178] D.C. Wdcox. A half a century historical review of the kw model. AIAA Paper 910615, AIAA 29th Aerospace Sciences Meeting, Reno, Nv,January 1991. [1791 F. Woodward. An improved method for the aercdynamic analysis of wmgbodytail configurations in subsonic and su rsonic flow, Part 1  the0 and application. N&A CR 2228 Pt.1, May 197? [MO] P. Woodward and P. Colella. ?e numerical simulation of twodimensional flutd flow with strong shocks. J. Comp. Phys., 54115173,1984. [181] K. Xu.L. Martinelli, and A. Jameson. Gaskinetic finite volume methods, fluxvector splitting and artiftcial diffusion. J. Comp. Phys., 1204865,1995.
f..Basic theory.
[ 1821 V. Yakhot and S.A. Orsza Renormalization
analysis of turbulence. Comp., 1:351, 1986.
oup
~ S C L
[1831 H.C. Yee. eaic and upwind TVD schemes. In P%c. th GAMM Conference on Numrical Methods in Fluid Mechanics. Gottinaen. September 1985.
SF
[ 1841 S. Yoon and A. Jameson. LowerU per Symmetric
GaussSeidel method for the Eu!er and NavierStokes equations. AIAA Paper 870600, AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.
I
Figure 21: NavierStokes Predictions for the F18 WingFuselage at Large Incident,
22b: NACA0012 Airfoil
22a: RAE2822 Airfoil
Figure 22: 0Topology Meshes, 160x32
'1  4
 3
 5
a 8
I  5
%:',
&
t&i
d*
a&
Figure 23: RAF,2822 Airfoil at Mach 0.750 and a=3.O0HCUSPScheme.
A
ra
Y
132
'1
P
wob
2 k C, after 35 Cycles.
24b Convergence.
Ci=0.3654, Cp0.0232.
Figure 2 4 NACA0012Airfoil at Mach 0.800 and C Z = ~ . ~ ~ ~ H Scheme. CUSP
4
f
B % % I
a
3

I:
I
25a: C, after 35 Cycles. Ci=0.3861, Cd=0.0582.
$  Wod
25b Convergence.
Figure 25: NACA0012Airfoil at Mach 0.850 and a=l.OoHCLJSPScheme.
133
a
*
3
e 4
/*+ '

... *..
;
2 1
si
**.
26b: 31.25% Span.
26a: 12.50% Span.
Cp0.2933, Cd=0.0274.
:/ *
.............. *+:*a.
c
3
:4*..
4:
I .
9
................:
*/
:
Ci=0.3139,Cd=0.0159.
.....
*.
t
26c: 50.00% Span.
Ci~0.3262,Cd=0.0089.
26d: 68.75% Span. Cl=0.3195,Cd=0.0026.
Figure 2 6 Onera M6 Wmg. Mach 0.840, Angle of Attack 3.06", 192x32~48Mesh. C~=0.3041,C~=0.0131 BCUSP scheme.
134
27a: Initial Wing CL=O.SOOI,C~=0.0210, a= 1.672"
27b 40 Design Iterations C~=0.5000, c&.Ol12, a=0.283'
Figure 27: Swept Wing Design Case (1). M=0.85, Fixed Lift Mode.Drag Reduction at CL=.^.
135
(1' UPPER SURFACE PRESSURE
UPPER SURFACE PRESSURE
28a: Initial Wing C~=0.5001,Cp0.0210,a= 1.672"
28b: 40 Design Iterations C~=0.5000,C~=0.0112,u=0.283'
Figure 28: Swept Wing Design Case (I), M=0.85, Fixed Lift Mode.Drag Reduction at C~=0.50
136
9
$11
5
a
2
.*
4
% 9 8.
'i
1 ,
'1 i 3 
29a: span station z=O.OO
29b span station 24.312
t
2%: span station ~ 4 . 6 2 5
29d: span station ~ 0 . 9 3 7
Figure 29: KO67 solution for initial wing.M=0.85,C~=0.4997, C0=0.G207, (r=1.970".
137
B
? 8 a
 
U"? 8
P
I.'
t f
f!
'1 !J
30a: span station z=O.OO
30b: span station z=0.312
8
$1
71B
U"?
'i ,j
 f!J
30c: span station z=0.625
30d: span station z=0.937
Figure 30: FL067 check on redesigned wing.M=0.85, C~=0.4992,C~=0.0094,a=0.300"
138
 e
s""*x++++++++++++++++++++++
;F
9
...............
.
+
8 
.""
9 
+
+:
+. t .
i
. t
c
t t
+
i
s
"
,"
. . 'I".
31b: span station z=0.312
 
I E
'1
++++
ii
31a: span station z=O.OO
SJ
+c .......... ".
a 
31c: span station z=0.625
'1 S i
31d: span station z=0.937
Figure 3 1: FL067 check on redesigned wing at a higher Mach number. M=0.86,C~=0.4988,C~=0.0097, (r=O.44O0.