tion of panel methods which could solve the linear flow models for arbitrarily ... supercritical airfoils [15]. Computational Fluid Dynamics (CFD) has...

5 downloads 4 Views 1MB Size

1

SUMMARY

der of approximations (p) has been successfully exploited both separately and in combination in the hp method [126]. A continuing obstacle to the treatment of conﬁgurations with complex geometry has been the problem of mesh generation. Several general techniques have been developed, including algebraic transformations and methods based on the solution of elliptic and hyperbolic equations. In the last few years methods using unstructured meshes have also begun to gain more general acceptance. The Dassault-INRIA group led the way in developing a ﬁnite element method for transonic potential ﬂow. They obtained a solution for a complete Falcon 50 as early as 1982 [25]. Euler methods for unstructured meshes have been the subject of intensive development by several groups since 1985 [110, 82, 81, 163, 14], and Navier-Stokes methods on unstructured meshes have also been demonstrated [117, 118, 11].

This paper presents a perspective on computational ﬂuid dynamics as a tool for aircraft design. It addresses the requirements for eﬀective industrial use, and trade-oﬀs between modelling accuracy and computational costs. Issues in algorithm design are discussed in detail, together with a uniﬁed approach to the design of shock capturing algorithms. Finally, the paper discusses the use of techniques drawn from control theory to determine optimal aerodynamic shapes. In the future multidisciplinary analysis and optimization should be combined to provide an integrated numerical design environment.

2

INTRODUCTION

Computational methods ﬁrst began to have a significant impact on aerodynamics analysis and design in the period of 1965-75. This decade saw the introduction of panel methods which could solve the linear ﬂow models for arbitrarily complex geometry in both subsonic and supersonic ﬂow [58, 147, 178]. It also saw the appearance of the ﬁrst satisfactory methods for treating the nonlinear equations of transonic ﬂow [123, 122, 63, 64, 43, 54], and the development of the hodograph method for the design of shock free supercritical airfoils [15].

Despite the advances that have been made, CFD is still not being exploited as eﬀectively as one would like in the design process. This is partly due to the long set-up and high costs, both human and computational of complex ﬂow simulations. The essential requirements for industrial use are: 1. assured accuracy 2. acceptable computational and human costs 3. fast turn around.

Computational Fluid Dynamics (CFD) has now matured to the point at which it is widely accepted as a key tool for aerodynamic design. Algorithms have been the subject of intensive development for the past two decades. The principles underlying the design and implementation of robust schemes which can accurately resolve shock waves and contact discontinuities in compressible ﬂows are now quite well established. It is also quite well understood how to design high order schemes for viscous ﬂow, including compact schemes and spectral methods. Adaptive reﬁnement of the mesh interval (h) and the or-

Improvements are still needed in all three areas. In particular, the ﬁdelity of modelling of high Reynolds number viscous ﬂows continues to be limited by computational costs. Consequently accurate and cost effective simulation of viscous ﬂow at Reynolds numbers associated with full scale ﬂight, such as the prediction of high lift devices, remains a challenge. Several routes are available toward the reduction of computational costs, including the reduction of mesh requirements by the use of higher order schemes, im1

proved convergence to a steady state by sophisticated acceleration methods, fast inversion methods for implicit schemes, and the exploitation of massively parallel computers.

dissipated by viscosity. The ratio of the length scale of the global ﬂow to that of the smallest persisting 3 eddies is of the order Re 4 , where Re is the Reynolds number, 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 9 of Re 4 cells would be required. This is beyond the range of any current or foreseeable computer. Consequently mathematical models with varying degrees of simpliﬁcation have to be introduced in order to make computational simulation of ﬂow feasible and produce viable and cost-eﬀective methods.

Another factor limiting the eﬀective use of CFD is the lack of good interfaces to computer aided design (CAD) systems. The geometry models provided by existing CAD systems often fail to meet the requirements of continuity and smoothness needed for ﬂow simulation, with the consequence that they must be modiﬁed before they can be used to provide the input for mesh generation. This bottleneck, which impedes the automation of the mesh generation process, needs to be eliminated, and the CFD software should be fully integrated in a numerical design environment. In addition to more accurate and cost-eﬀective ﬂow prediction methods, better optimizations methods are also needed, so that not only can designs be rapidly evaluated, but directions of improvement can be identiﬁed. Possession of techniques which result in a faster design cycle gives a crucial advantage in a competitive environment.

Figure 1 (supplied by Pradeep Raj) indicates a hierarchy of models at diﬀerent levels of simpliﬁcation which have proved useful in practice. Eﬃcient ﬂight is generally achieved by the use of smooth and streamlined shapes which avoid ﬂow separation and minimize viscous eﬀects, with the consequence that useful predictions can be made using inviscid models. Inviscid calculations with boundary layer corrections can provide quite accurate predictions of lift and drag when the ﬂow remains attached, but iteration between the inviscid outer solution and the inner boundary layer solution becomes increasingly diﬃcult with the onset of separation. Procedures for solving the full viscous equations are likely to be needed for the simulation of arbitrary complex separated ﬂows, which may occur at high angles of attack or with bluﬀ bodies. In order to treat ﬂows at high Reynolds numbers, one is generally forced to estimate turbulent eﬀects by Reynolds averaging of the ﬂuctuating components. This requires the introduction of a turbulence model. As the available computing power increases one may also aspire to large eddy simulation (LES) in which the larger scale eddies are directly calculated, while the inﬂuence of turbulence at scales smaller than the mesh interval is represented by a subgrid scale model.

A critical issue, examined in the next section, is the choice of mathematical models. What level of complexity is needed to provide suﬃcient accuracy for aerodynamic design, and what is the impact on cost and turn-around time? Section 3 addresses the design of numerical algorithms for ﬂow simulation. Section 4 presents the results of some numerical calculations which require moderate computer resources and could be completed with the fast turn-around required by industrial users. Section 5 discusses automatic design procedures which can be used to produce optimum aerodynamic designs. Finally, Section 7 oﬀers an outlook for the future.

3

3.1

THE COMPLEXITY OF FLUID FLOW AND MATHEMATICAL MODELLING The Hierarchy of Mathematical Models

Many critical phenomena of ﬂuid ﬂow, such as shock waves and turbulence, are essentially nonlinear. They also exhibit extreme disparities of scales. While the actual thickness of a shock wave is of the order of a mean free path of the gas particles, on a macroscopic scale its thickness is essentially zero. In turbulent ﬂow energy is transferred from large scale motions to progressively smaller eddies until the scale becomes so small that the motion is

Figure 1: Hierarchy of Fluid Flow Models

2

3.2

Computational Costs

Computational costs vary drastically with the choice of mathematical model. Panel methods can be effectively used to solve the linear potential ﬂow equation with higher-end personal computers (with an Intel 80486 microprocessor, for example). Studies of the dependency of the result on mesh reﬁnement, performed by this author and others, have demonstrated that inviscid transonic potential ﬂow 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 enough to obtain a converged result. Consequently airfoil calculations can be performed in seconds on a Cray YMP, and can also be performed on 486-class personal computers. Correspondingly accurate three-dimensional inviscid calculations can be performed for a wing on a mesh, say with 192×32×48 = 294, 912 cells, in about 5 minutes on a single processor Cray YMP, or less than a minute with eight processors, or in 1 or 2 hours on a workstation such as a Hewlett Packard 735 or an IBM 560 model.

32 cells 32 cells in the boundary layer

512 cells around the wing to limit the mesh aspect ratio (to about 1000)

Surface Mesh

256 cells spanwise

Total: 512 x 64 x 256 = 8 388 608 cells

Figure 2: Mesh Requirements for a Viscous Simulation

pose that a conservative estimate of the size of eddies in a boundary layer that ought to be resolved is 1/5 of the boundary layer thickness. Assuming that 10 points are needed to resolve 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 for a wing with an aspect ratio of 10, 50,000 intervals will be needed in the spanwise direction. Thus, of the order of 50×5, 000×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 total time equal to the time required for waves to travel three times the length of the chord, of the order of 15,000 time steps would be needed. Performance beyond the teraﬂop (1012 operations per second) will be needed to attempt calculations of this nature, which also have an information content far beyond what is needed for enginering analysis and design. The designer does not need to know the details of the eddies in the boundary layer. The primary purpose of such calculations is to improve the calculation of averaged quantities such as skin friction, and the prediction of global behavior such as the onset of separation. The main current use of Navier-Stokes and large eddy simulations is to try to gain an improved insight into the physics of turbulent ﬂow, which may in turn lead to the development of more comprehensive and reliable

Viscous simulations at high Reynolds numbers require vastly greater resources. Careful twodimensional studies of mesh requirements have been carried out at Princeton by Martinelli [114]. He found that on the order of 32 mesh intervals were needed to resolve a turbulent boundary layer, in addition to 32 intervals between the boundary layer and the far ﬁeld, leading to a total of 64 intervals. In order to prevent degradations in accuracy and convergence due to excessively large aspect ratios (in excess of 1,000) in the surface mesh cells, the chordwise resolution must also be increased to 512 intervals. Reasonably accurate solutions can be obtained in a 512×64 mesh in 100 multigrid cycles. Translated to three dimensions, this would imply the need for meshes with 5–10 million cells (for example, 512×64×256 = 8,388,608 cells as shown in Figure 2). When simulations are performed on less ﬁne meshes with, say, 500,000 to 1 million cells, it is very hard to avoid mesh dependency in the solutions as well as sensitivity to the turbulence model. A typical algorithm requires of the order of 5,000 ﬂoating point operations per mesh point in one multigrid iteration. With 10 million mesh points, the operation count is of the order of 0.5×1011 per cycle. Given a computer capable of sustaining 1011 operations per second (100 gigaﬂops), 200 cycles could then be performed in 100 seconds. Simulations of unsteady viscous ﬂows (ﬂutter, buﬀet) would be likely to require 1,000–10,000 time steps. A further progression to large eddy simulation of complex conﬁgurations would require even greater resources. The following estimate is due to W.H. Jou [90]. Sup3

turbulence models.

3.3

timately rests with industry. Aircraft and spacecraft designs normally pass through the three phases of conceptual design, preliminary design, and detailed design. Correspondingly, the appropriate CFD models will vary in complexity. In the conceptual and preliminary design phases, the emphasis will be on relatively simple models which can give results with very rapid turn-around and low computer costs, in order to evaluate alternative conﬁgurations and perform quick parametric studies. The detailed design stage requires the most complete simulation that can be achieved with acceptable cost. In the past, the low level of conﬁdence that could be placed on numerical predictions has forced the extensive use of wind tunnel testing at an early stage of the design. This practice was very expensive. The limited number of models that could be fabricated also limited the range of design variations that could be evaluated. It can be anticipated that in the future, the role of wind tunnel testing in the design process will be more one of veriﬁcation. Experimental research to improve our understanding of the physics of complex ﬂows will continue, however, to play a vital role.

Turbulence Modelling

It is doubtful whether a universally valid turbulence model, capable of describing all complex ﬂows, could be devised [52]. Algebraic models [30, 9] have proved fairly satisfactory for the calculation of attached and slightly separated wing ﬂows. These models rely on the boundary layer concept, usually incorporating separate formulas for the inner and outer layers, and they require an estimate of a length scale which 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 in the Baldwin-Lomax model, can lead to ambiguities in internal ﬂows, and also in complex vortical ﬂows over slender bodies and highly swept or delta wings [40, 115]. The Johnson-King model [88], which allows for non-equilibrium eﬀects through the introduction of an ordinary diﬀerential equation for the maximum shear stress, has improved the prediction of ﬂows with shock induced separation [148, 91]. Closure models depending on the solution of transport equations are widely accepted for industrial applications. These models eliminate the need to estimate a length scale by detecting the edge of the boundary layer. Eddy viscosity models typically use two equations for the turbulent kinetic energy k and the dissipation rate , or a pair of equivalent quantities [89, 177, 160, 1, 121, 35]. Models of this type generally tend to present diﬃculties in the region very close to the wall. They also tend to be badly conditioned for numerical solution. The k − l model [154] is designed to alleviate this problem by taking advantage of the linear behaviour of the length scale l near the wall. In an alternative approach to the design of models which are more amenable to numerical solution, new models requiring the solution of one transport equation have recently been introduced [10, 159]. The performance of the algebraic models remains competitive for wing ﬂows, but the one- and two-equation models show promise for broader classes of ﬂows. In order to achieve greater universality, research is also being pursued on more complex Reynolds stress transport models, which require the solution of a larger number of transport equations.

4 4.1

CFD ALGORITHMS Diﬃculties of Flow Simulation

The computational simulation of ﬂuid ﬂow presents a number of severe challenges for algorithm design. At the level of inviscid modeling, the inherent nonlinearity of the ﬂuid ﬂow equations leads to the formation of singularities such as shock waves and contact discontinuities. Moreover, the geometric conﬁgurations of interest are extremely complex, and generally contain sharp edges which lead to the shedding of vortex sheets. Extreme gradients near stagnation points or wing tips may also lead to numerical errors that can have global inﬂuence. Numerically generated entropy may be convected from the leading edge for example, causing the formation of a numerically induced boundary layer which can lead to separation. The need to treat exterior domains of inﬁnite extent is also a source of diﬃculty. Boundary conditions imposed at artiﬁcial outer boundaries may cause reﬂected waves which signiﬁcantly interfere with the ﬂow. When viscous eﬀects are also included in the simulation, the extreme diﬀerence of the scales in the viscous boundary layer and the outer ﬂow, which is essentially inviscid, is another source of diﬃculty, forcing the use of meshes with extreme variations in mesh interval. For these reasons CFD, has been a driving force for the development of numerical algorithms.

Another direction of research is the attempt to devise more rational models via renormalization group (RNG) theory [181, 155]. Both algebraic and twoequation k − models devised by this approach have shown promising results [116]. The selection of suﬃciently accurate mathematical models and a judgment of their cost eﬀectiveness ul4

4.2

Structured Meshes

and

Unstructured

the beginning of the century by Voronoi [174], and the moving front method [111].

4.3

The algorithm designer faces a number of critical decisions. The ﬁrst choice that must be made is the nature of the mesh used to divide the ﬂow ﬁeld into discrete subdomains. The discretization procedure must allow for the treatment of complex conﬁgurations. The principal alternatives are Cartesian meshes, body-ﬁtted curvilinear meshes, and unstructured tetrahedral meshes. Each of these approaches has advantages which have led to their use. The Cartesian mesh minimizes the complexity of the algorithm at interior points and facilitates the use of high order discretization procedures, at the expense of greater complexity, and possibly a loss of accuracy, in the treatment of boundary conditions at curved surfaces. This diﬃculty may be alleviated by using mesh reﬁnement procedures near the surface. With their aid, schemes which use Cartesian meshes have recently been developed to treat very complex conﬁgurations [120, 149, 22, 94].

Finite Diﬀerence, Finite Volume, and Finite Element Schemes

Associated with choice of mesh type is the formulation of the discretization procedure for the equations of ﬂuid ﬂow, which can be expressed as diﬀerential conservation laws. In the Cartesian tensor notation, let xi be the coordinates, p, ρ, T , and E the pressure, density, temperature, and total energy, and ui the velocity components. Using the convention that summation over j = 1, 2, 3 is implied by a repeated subscript j, each conservation equation has the form ∂fj ∂w + = 0. ∂t ∂xj

(1)

For the mass equation w = ρ, fj = ρuj . For the i momentum equation

Body-ﬁtted meshes have been widely used and are particularly well suited to the treatment of viscous ﬂow because they readily allow the mesh to be compressed near the body surface. With this approach, the problem of mesh generation itself has proved to be a major pacing item. The most commonly used procedures are algebraic transformations [7, 44, 46, 156], methods based on the solution of elliptic equations, pioneered by Thompson [169, 170, 157, 158], and methods based on the solution of hyperbolic equations marching out from the body [161]. In order to treat very complex conﬁgurations it generally proves expedient to use a multiblock [176, 150] procedure, with separately generated meshes in each block, which may then be patched at block faces, or allowed to overlap, as in the Chimera scheme [19, 20]. While a number of interactive software systems for grid generation have been developed, such as EAGLE, GRIDGEN, and ICEM, the generation of a satisfactory grid for a very complex conﬁguration may require months of eﬀort.

wi = ρui , fij = ρui uj + pδij − σij , where σij is the viscous stress tensor. For the energy equation w = ρE, fj = (ρE + p)uj − σjk uk − κ

∂T , ∂xj

where κ is the coeﬃcient of heat conduction. The pressure is related to the density and energy by the equation of state 1 p = (γ − 1)ρ E − ui ui (2) 2 in which γ is the ratio of speciﬁc heats. In the Navier-Stokes equations the viscous stresses are assumed to be linearly proportional to the rate of strain, or ∂ui ∂uk ∂uj σij = µ + + λδij , (3) ∂xj ∂xi ∂xk where µ and λ are the coeﬃcients of viscosity and bulk viscosity, and usually λ = −2µ/3.

The alternative is to use an unstructured mesh in which the domain is subdivided into tetrahedra. This in turn requires the development of solution algorithms capable of yielding the required accuracy on unstructured meshes. This approach has been gaining acceptance, as it is becoming apparent that it can lead to a speed-up and reduction in the cost of mesh generation that more than oﬀsets the increased complexity and cost of the ﬂow simulations. Two competing procedures for generating triangulations which have both proved successful are Delaunay triangulation [41, 11], based on concepts introduced at

The ﬁnite diﬀerence method, which requires the use of a Cartesian or a structured curvilinear mesh, directly approximates the diﬀerential operators appearing in these equations. In the ﬁnite volume method [112], the discretization is accomplished by dividing the domain of the ﬂow into a large number of small subdomains, and applying the conservation laws in the integral form ∂ wdV + f · dS = 0. ∂t Ω ∂Ω 5

Here F is the ﬂux appearing in equation (1) and dS is the directed surface element of the boundary ∂Ω of the domain Ω. The use of the integral form has the advantage that no assumption of the diﬀerentiability of the solutions is implied, with the result that it remains a valid statement for a subdomain containing a shock wave. In general the subdomains could be arbitrary, but it is convenient to use either hexahedral cells in a body conforming curvilinear mesh or tetrahedrons in an unstructured mesh.

(4) has the advantages that it is valid for both structured and unstructured meshes, and that it assures that a uniform ﬂow exactly satisﬁes the equations, because faces S = 0 for a closed control volume. Finite diﬀerence schemes do not necessarily satisfy this constraint because of the discretization errors ∂ξi and the inversion of the transforin evaluating ∂x j mation matrix. A cell-vertex ﬁnite volume scheme can be derived by taking the union of the cells surrounding a given vertex as the control volume for that vertex [55, 71, 139]. In equation (4), V is now the sum of the volumes of the surrounding cells, while the ﬂux balance is evaluated over the outer faces of the polyhedral control volume. In the absence of upwind biasing the ﬂux vector is evaluated by averaging 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 ﬁnite element method. Whereas the ﬁnite diﬀerence and ﬁnite volume methods approximate the diﬀerential and integral operators, the ﬁnite element method proceeds by inserting an approximate solution into the exact equations. On multiplying by a test function φ and integrating by parts over space, one obtains the weak form

Alternative discretization schemes may be obtained by storing ﬂow variables at either the cell centers or the vertices. These variations are illustrated in Figure 3 for the two-dimensional case. With a cellcentered scheme the discrete conservation law takes the form

o

o

3a: Cell Centered Scheme.

∂ ∂t

o

φwdΩ =

Ω

f ·∇φdΩ −

Ω

φf ·dS (6) ∂Ω

which is also valid in the presence of discontinuities in the ﬂow. In the Galerkin method the approximate solution is expanded in terms of the same family of functions as those from which the test functions are drawn. By choosing test functions with local support, separate equations are obtained for each node. For example, if a tetrahedral mesh is used, and φ is piecewise linear, with a nonzero value only at a single node, the equations at each node have a stencil which contains only the nearest neighbors. In this case the ﬁnite element approximation corresponds closely to a ﬁnite volume scheme. If a piecewise linear approximation to the ﬂux f is used in the evaluation of the integrals on the right hand side of equation (6), these integrals reduce to formulas which are identical to the ﬂux balance of the ﬁnite volume scheme.

o

3b: Vertex Scheme. Figure 3: Structured and Unstructured Discretizations. d f ·S = 0, (4) wV + dt faces where V is the cell volume, and f is now a numerical estimate of the ﬂux vector through each face. f may be evaluated from values of the ﬂow variables in the cells separated by each face, using upwind biasing to allow for the directions of wave propagation. With hexahedral cells, equation (4) is very similar to a ﬁnite diﬀerence scheme in curvilinear coordinates. Under a transformation to curvilinear coordinates ξj , equation (1) becomes ∂ ∂ ∂ξi (Jw) + fj = 0, J (5) ∂t ∂ξi ∂xj

Thus the ﬁnite diﬀerence and ﬁnite volume methods lead to essentially similar schemes on structured meshes, while the ﬁnite volume method is essentially equivalent to a ﬁnite element method with linear elements when a tetrahedral mesh is used. Provided that the ﬂow equations are expressed in the conservation law form (1), all three methods lead to an exact cancellation of the ﬂuxes through interior cell boundaries, so that the conservative property of the equations is preserved. The important role of this

where J is the Jacobian determinant of the transfor∂ξi i mation matrix ∂x ∂ξj . The transformed ﬂux J ∂xj fj corresponds to the dot product of the ﬂux f with a ∂ξi , while J represents the transvector face area J ∂x j formation of the cell volume. The ﬁnite volume form 6

property in ensuring correct shock jump conditions was pointed out by Lax and Wendroﬀ [97].

4.4 4.4.1

of a solution of this equation does not increase, provided that any discontinuity appearing in the solution satisﬁes an entropy condition [96]. Harten proposed that diﬀerence schemes ought to be designed so that the discrete total variation cannot increase [56]. If the end values are ﬁxed, the total variation can be expressed as

T V (u) = 2 maxima − minima .

Non-oscillatory Shock Capturing Schemes Local Extremum Diminishing (LED) Schemes

Thus a LED scheme is also total variation diminishing (TVD). Positivity conditions of the type expressed in equations (7) and (8) lead to diagonally dominant schemes, and are the key to the elimination of improper oscillations. The positivity conditions may be realized by the introduction of diﬀusive terms or by the use of upwind biasing in the discrete scheme. Unfortunately, they may also lead to severe restrictions on accuracy unless the coeﬃcients have a complex nonlinear dependence on the solution.

The discretization procedures which have been described in the last section lead to nondissipative approximations to the Euler equations. Dissipative terms may be needed for two reasons. The ﬁrst is the possibility of undamped oscillatory modes. The second reason is the 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 uniﬁed approach to the construction of nonoscillatory schemes via the introduction of controlled diﬀusive and antidiﬀusive terms. This is the line adhered to in the author’s own work.

4.4.2

Following the pioneering work of Godunov [51], 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, 165, 5, 68, 182, 62, 179, 13, 12, 11]. If the one-dimensional scalar conservation law ∂ ∂v + f (v) = 0 (9) ∂t ∂x is represented by a three point scheme

The development of non-oscillatory schemes has been a prime focus of algorithm research for compressible ﬂow. Consider a general semi-discrete scheme of the form d vj = cjk (vk − vj ). dt

(7)

k=j

A maximum cannot increase and a minimum cannot decrease if the coeﬃcients cjk are non-negative, since at a maximum vk − vj ≤ 0, and at a minimum vk − vj ≥ 0. Thus the condition cjk ≥ 0, k = j

Artiﬁcial Diﬀusion and Upwinding

dvj = c+ (vj+1 − vj ) + c− (vj−1 − vj ) , j+ 12 j− 12 dt the scheme is LED if

(8)

≥ 0, c− ≥ 0. c+ j+ 1 j− 1 2

is suﬃcient to ensure stability in the maximum norm. Moreover, if the scheme has a compact stencil, so that cjk = 0 when j and k are not nearest neighbors, a local maximum cannot increase and local minimum cannot decrease. This local extremum diminishing (LED) property prevents the birth and growth of oscillations. The one-dimensional conservation law ∂ ∂u + f (u) = 0 ∂t ∂x

2

(10)

A conservative semidiscrete approximation to the one-dimensional conservation law can be derived by subdividing the line into cells. Then the evolution of the value vj in the jth cell is given by ∆x

dvj + hj+ 12 − hj− 12 = 0, dt

(11)

where hj+ 12 is an estimate of the ﬂux between cells j and j + 1. The simplest estimate is the arithmetic average (fj+1 + 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

provides a useful model for analysis. In this case waves are propagated with a speed a(u) = ∂f ∂u , and the solution is constant along the characteristics dx dt = a(u). Thus the LED property is satisﬁed. In fact the total variation ∞ ∂u dx T V (u) = −∞ ∂x

hj+ 12 =

1 (fj+1 + fj ) − αj+ 12 (vj+1 − vj ) . 2

(12)

In order to estimate the required value of the coeﬃcient αj+ 12 , let aj+ 12 be a numerical estimate of the 7

wave speed

∂f ∂u ,

aj+ 12 =

fj+1 −fj vj+1 −vj ∂f ∂v v=vj

Then hj+ 12 − hj− 12

if vj+1 = vj if vj+1 = vj

.

accurate in the neighborhood of an extremum. It proves useful in the following development to introduce the concept of essentially local extremum diminishing (ELED) schemes. These are deﬁned to be schemes which satisfy the condition that in the limit as the mesh width ∆x → 0, local maxima are nonincreasing, and local minima are non-decreasing.

(13)

=

1 aj+ 12 − αj+ 12 ∆vj+ 12 2 1 aj− 12 + αj− 12 ∆vj− 12 , + 2 +

4.4.3

where

Higher order non-oscillatory schemes can be derived by introducing anti-diﬀusive terms in a controlled manner. An early attempt to produce a high resolution scheme by this approach is the JamesonSchmidt-Turkel (JST) scheme [85]. Suppose that anti-diﬀusive terms are introduced by subtracting neighboring diﬀerences to produce a third order diffusive ﬂux

∆vj+ 12 = vj+1 − vj , and the LED condition (10) is satisﬁed if 1 αj+ 12 ≥ aj+ 12 . 2

High Resolution Switched Schemes: Jameson-SchmidtTurkel (JST) Scheme

(14)

If one takes

1 aj+ 12 , 2 one obtains the ﬁrst order upwind scheme fj if aj+ 12 > 0 hj+ 12 = . fj+1 if aj+ 12 < 0 αj+ 12 =

dj+ 1 = αj+ 1 2

2

∆vj+ 1 − 2

1 ∆vj+ 3 + ∆vj− 1 2 2 2

, (15)

3

∂ which is an approximation to 12 α∆x3 ∂x 3 . The positivity condition (8) is violated by this scheme. It proves that it generates substantial oscillations in the vicinity of shock waves, which can be eliminated by switching locally to the ﬁrst order scheme. The JST scheme therefore introduces blended diﬀusion of the form

This is the least diﬀusive ﬁrst order scheme which satisﬁes 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 potential ﬂow also involved the use of upwind biasing. This was ﬁrst introduced by Murman and Cole to treat the transonic small disturbance equation [123].

=

dj+ 12

−

(2)

+j+ 1 ∆vj+ 12 (16) 2

(4) j+ 1 ∆vj+ 32 − 2∆vj+ 12 + ∆vj− 12 , 2

(2)

Another important requirement of discrete schemes is that they should exclude 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, 168, 128, 131]. In the case that the numerical ﬂux function is strictly convex, Aiso has recently proved [2] that it is suﬃcient that 1 αj+ 12 > max aj+ 12 , sign(vj+1 − vj ) 2

(4)

The idea is to use variable coeﬃcients j+ 1 and j+ 1 2 2 which produce a low level of diﬀusion in regions where the solution is smooth, but prevent oscilla(2) tions near discontinuities. If j+ 1 is constructed so 2

that it is of order ∆x2 where the solution is smooth, (4) while j+ 1 is of order unity, both terms in dj+ 12 will 2

be of order ∆x3 .

The JST scheme has proved very eﬀective in practice in numerous calculations of complex steady ﬂows, and conditions under which it could be a total variation diminishing (TVD) scheme have been examined by Swanson and Turkel [164]. An alternative statement of suﬃcient conditions on the coeﬃcients (2) (4) j+ 1 and j+ 1 for the JST scheme to be LED is as 2 2 follows:

for > 0. Thus the numerical viscosity should be rounded out and not allowed to reach zero at a point where the wave speed a(u) = ∂f ∂u approaches zero. This justiﬁes, for example, Harten’s entropy ﬁx [56]. Higher order schemes can be constructed by introducing higher order diﬀusive terms. Unfortunately these have larger stencils and coeﬃcients of varying sign which are not compatible with the conditions (8) for a LED scheme, and it is known that schemes which satisfy these conditions are at best ﬁrst order

Theorem 1 (Positivity of the JST scheme) Suppose that whenever either vj+1 or vj is an extremum the coeﬃcients of the JST scheme satisfy 1 (4) (2) j+ 1 ≥ αj+ 12 , j+ 1 = 0. (17) 2 2 2 8

Then the JST scheme is local extremum diminishing (LED).

P4. L(u, v) = 0 if u and v have opposite signs: otherwise L(u, v) has the same sign as u and v.

Proof: We need only consider the rate of change of v at extremal points. Suppose that vj is an extremum. Then (4)

Properties (P1–P3) are natural properties of an average. Property (P4) is needed for the construction of a LED or TVD scheme.

(4)

j+ 1 = j− 1 = 0, 2

2

It is convenient to introduce the notation

and the semi-discrete scheme (11) reduces to 1 dvj (2) = j+ 1 − aj+ 12 ∆vj+ 12 ∆x 2 dt 2 1 (2) − j− 1 + aj− 12 ∆vj− 12 , 2 2

φ(r) = L(1, r) = L(r, 1), where according to (P4) φ(r) ≥ 0. It follows from (P2) on setting α = u1 or v1 that u

v

u=φ v. L(u, v) = φ u v

and each coeﬃcient has the required sign. ✷ (2) j− 1 2

Also it follows on setting v = 1 and u = r that 1 φ(r) = rφ . r

(4) j− 1 2

In order to construct and with the desired properties deﬁne u−v q 0 or v = 0 |u|+|v| if u = R(u, v) = (18) 0 if u = v = 0,

Thus, if there exists r < 0 for which φ(r) > 0, then φ 1r < 0. The only way to ensure that φ(r) ≥ 0 is to require φ(r) = 0 for all r < 0, corresponding to property (P4).

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

Now one deﬁnes the diﬀusive ﬂux for a scalar conservation law as

Qj = R(∆vj+ 12 , ∆vj− 12 ), Qj+ 12 = max(Qj , Qj+1 ).

dj+ 1 = αj+ 1 2

and

∆vj+ 1 − L ∆vj+ 3 , ∆vj− 1

2

2

Set

(2)

(4)

j+ 1 = αj+ 1 Qj+ 1 , j+ 1 2

2

2

2

1 = αj+ 1 (1 − Qj+ 1 ). (19) 2 2 2

r+ =

2

∆vj+ 32 ∆vj− 12

, r− =

∆vj− 32 ∆vj+ 12

.

2

(20)

.

and 4.4.4

L(∆vj+ 32 , ∆vj− 12 ) =

Symmetric Limited Positive (SLIP) Scheme

L(∆vj− 32 , ∆vj+ 12 ) = Then,

An alternative route to high resolution without oscillation is to introduce ﬂux limiters to guarantee the satisfaction of the positivity condition (8). The use of limiters dates back to the work of Boris and Book [23]. A particularly simple way to introduce limiters, proposed by the author in 1984 [68], is to use ﬂux limited dissipation. In this scheme the third order diﬀusion deﬁned by equation (15) is modiﬁed by the insertion of limiters which produce an equivalent three point scheme with positive coeﬃcients. The original scheme [68] can be improved in the following manner so that less restrictive ﬂux limiters are required. Let L(u, v) be a limited average of u and v with the following properties:

∆x

dvj dt

−

=

αj+ 1 − 2

αj− 1 + 2

φ(r+ )∆vj− 12

φ(r− )∆vj+ 12 .

1 a 1 + αj− 1 φ(r − ) ∆vj+ 1 2 2 2 j+ 2

1 a 1 + αj+ 1 φ(r + ) ∆vj− 1 (21) . 2 2 2 j− 2

Thus the scheme satisﬁes the LED condition if αj+ 12 ≥ 12 aj+ 12 for all j, and φ(r) ≥ 0, which is assured by property (P4) on L. At the same time it follows from property (P3) that the ﬁrst order diffusive ﬂux is canceled when ∆v is smoothly varying and of constant sign. Schemes constructed by this formulation will be referred to as symmetric limited positive (SLIP) schemes. This result may be summarized as

P1. L(u, v) = L(v, u)

Theorem 2 (Positivity of the SLIP scheme) Suppose that the discrete conservation law (11) contains a limited diﬀusive ﬂux as deﬁned by equation (20). Then the positivity condition (14), together

P2. L(αu, αv) = αL(u, v) P3. L(u, u) = u 9

with the properties (P1–P4) for limited averages, are suﬃcient to ensure satisfaction of the LED principle that a local maximum cannot increase and a local minimum cannot decrease. ✷

4.4.5

The limiters deﬁned by the formula (22) have the disadvantage that they are active at a smooth extrema, reducing the local accuracy of the scheme to ﬁrst order. In order to prevent this, the SLIP scheme can be relaxed to give an essentially local extremum diminishing (ELED) scheme which is second order accurate at smooth extrema by the introduction of a threshold in the limited average. Therefore redeﬁne D(u, v) as q u−v , (23) D(u, v) = 1 − r max(|u| + |v| , ∆x )

A variety of limiters may be deﬁned which meet the requirements of properties (P1–P4). Deﬁne S(u, v) =

1 {sign(u) + sign(v)} 2

which vanishes is u and v have opposite signs. Then two limiters which are appropriate are the following well-known schemes:

where r = 32 , q ≥ 2. This reduces to the previous deﬁnition if |u| + |v| > ∆xr .

1. Minmod: L(u, v) = S(u, v) min(|u|, |v|)

In any region where the solution is smooth, ∆vj+ 32 − ∆vj− 12 is of order ∆x2 . In fact if there is a smooth extremum in the neighborhood of vj or vj+1 , a Taylor series expansion indicates that ∆vj+ 32 , ∆vj+ 12 and ∆vj− 12 are each individually of order ∆x2 , since dv dx = 0 at the extremum. It may be veriﬁed that second order accuracy is preserved at a smooth extremum if q ≥ 2. On hand the limiter acts the other in the usual way if ∆vj+ 32 or ∆vj− 32 > ∆xr , and it may also be veriﬁed that in the limit ∆x → 0 local maxima are non increasing and local minima are non decreasing [79]. Thus the scheme is essentially local extremum diminishing (ELED).

2. Van Leer: L(u, v) = S(u, v)

2|u||v| . |u| + |v|

In order to produce a family of limiters which contains these as special cases it is convenient to set L(u, v) =

1 D(u, v)(u + v), 2

where D(u, v) is a factor which should deﬂate the arithmetic average, and become zero if u and v have opposite signs. Take u − v q , (22) D(u, v) = 1 − R(u, v) = 1 − |u| + |v|

The eﬀect 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 updating scheme. In a scheme recently proposed by Venkatakrishnan a threshold is introduced precisely for this purpose [173].

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) = 0 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 equivalent to Van Leer’s limiter. By increasing q one can generate a sequence of limited averages which approach a limit deﬁned by the arithmetic mean truncated to zero when u and v have opposite signs.

4.4.6

When the terms are regrouped, it can be seen that with this limiter the SLIP scheme is exactly equivalent to the JST scheme, with the switch is deﬁned as

Qj+ 12 = R ∆vj+ 32 , ∆vj+ 12 (2)

j+ 1 2

(4)

j+ 1 2

Essentially Local Extremum Diminishing (ELED) Scheme with Soft Limiter

Upstream Limited Positive (USLIP) Schemes

By adding the anti-diﬀusive correction purely from the upstream side one may derive a family of upstream limited positive (USLIP) schemes. Corresponding to the original SLIP scheme deﬁned by equation (20), a USLIP scheme is obtained by setting

dj+ 12 = αj+ 12 ∆vj+ 12 − L ∆vj+ 12 , ∆vj− 12

= αj+ 12 Qj+ 12

= αj+ 12 1 − Qj+ 12 .

if aj+ 12 > 0, or

This formulation thus uniﬁes the JST and SLIP schemes.

dj+ 12 = αj+ 12 10

∆vj+ 12 − L ∆vj+ 12 , ∆vj+ 32

if aj+ 12 < 0. If αj+ 12 = 12 aj+ 12 one recovers a standard high resolution upwind scheme in semidiscrete form. Consider the case that aj+ 12 > 0 and aj− 12 > 0. If one sets r+ =

∆vj+ 12 ∆vj− 12

, r− =

the scheme reduces to

∆x

dvj 1 =− dt 2

∆vj− 32 ∆vj− 12

Following Roe’s derivation, let Aj+ 12 be a mean value Jacobian matrix exactly satisfying the condition

,

fj+1 − fj = Aj+ 12 (wj+1 − wj ).

φ(r + )aj+ 1 + 2 − φ(r − ) aj− 1 2

where now the ﬂux diﬀerence fj+1 − fj is split. The corresponding diﬀusive ﬂux is

1 + − . ∆fj+ 1 − ∆fj+ dj+ 12 = 1 2 2 2

2

Aj+ 12 may be calculated by substituting the weighted averages

∆vj− 1 . 2

√

To assure the correct sign to satisfy the LED criterion the ﬂux limiter must now satisfy the additional constraint that φ(r) ≤ 2.

u=

√ √ √ ρj+1 uj+1 + ρj uj ρj+1 Hj+1 + ρj Hj ,H = √ √ √ √ ρj+1 + ρj ρj+1 + ρj (27)

into the standard formulas for the Jacobian ma∂f trix A = ∂w . A splitting according to characteristic ﬁelds is now obtained by decomposing Aj+ 12 as

The USLIP formulation is essentially equivalent to standard upwind schemes [130, 165]. Both the SLIP and USLIP constructions can be implemented on unstructured meshes [75, 79]. The anti-diﬀusive terms are then calculated by taking the scalar product of the vectors deﬁning an edge with the gradient in the adjacent upstream and downstream cells. 4.4.7

Aj+ 12 = T ΛT −1,

Aj+ 12 = T |Λ| T −1

where Steger and Warming [162] ﬁrst showed how to generalize the concept of upwinding to the system of conservation laws

and |Λ| is the diagonal matrix containing the absolute values of the eigenvalues.

(24) 4.4.8

by the concept of ﬂux splitting. Suppose that the + ∂f − ﬂux is split as f = f + + f − where ∂f ∂w and ∂w have positive and negative eigenvalues. Then the ﬁrst order upwind scheme is produced by taking the numerical ﬂux to be

This can be expressed in viscosity form as =

=

1 + 1 + fj+1 + fj+ − fj+1 − fj+ 2 2 1 − 1 − − fj+1 − fj− + fj+1 + fj + 2 2 1 (fj+1 + fj ) − dj+ 12 , 2

dj+ 12 =

1 ∆(f + − f − )j+ 12 . 2

1 α 1 ∆wj+ 12 . 2 j+ 2

(29)

The JST scheme with scalar diﬀusive ﬂux captures shock waves with about 3 interior points, and it has been widely used for transonic ﬂow calculations because it is both robust and computationally inexpensive.

where the diﬀusive ﬂux is dj+ 12 =

Alternative Splittings

Characteristic splitting has the advantages that it introduces the minimum amount of diﬀusion to exclude the growth of local extrema of the characteristic variables, and that with the Roe linearization it allows a discrete shock structure with s single interior point. To reduce the computational complexity one may replace |A| by αI where if α is at least equal to the spectral radius max |λ(A)|, then the positivity conditions will still be satisﬁed. Then the ﬁrst order scheme simply has the scalar diﬀusive ﬂux

− . hj+ 12 = fj+ + fj+1

hj+ 12

(28)

where the columns of T are the eigenvectors of Aj+ 12 , and Λ is a diagonal matrix of the eigenvalues. Now the corresponding diﬀusive ﬂux is 1 Aj+ 12 (wj+1 − wj ), 2

Systems of Conservation Laws: Flux Splitting and Flux-Diﬀerence Splitting

∂ ∂w + f (w) = 0 ∂t ∂x

(26)

(25)

Roe derived the alternative formulation of ﬂux difference splitting [146] by distributing the corrections due to the ﬂux diﬀerence in each interval upwind and downwind to obtain

An intermediate class of schemes can be formulated by deﬁning the ﬁrst order diﬀusive ﬂux as a combination of diﬀerences of the state and ﬂux vectors dj+ 1 =

dwj + (fj+1 − fj )− + (fj − fj−1 )+ = 0, ∆x dt

2

11

1 ∗ 1 α 1 c (wj+1 − wj ) + βj+ 1 (fj+1 − fj ) , 2 2 j+ 2 2 (30)

where the factor c is included in the ﬁrst term to make α∗j+ 1 and βj+ 12 dimensionless. Schemes of this 2 class are fully upwind in supersonic ﬂow if one takes α∗j+ 1 = 0 and βj+ 12 = sign(M ) when the absolute 2 value of the Mach number M exceeds 1. The ﬂux vector f can be decomposed as f = uw + fp ,

where

wL

wL wA

j-2

(31)

j-1 j

wR

0 fp = p . up

(32)

j+1

wR j+2

Figure 4: Shock structure for single interior point.

Then fj+1 − fj = u ¯ (wj+1 − wj ) + w ¯ (uj+1 − uj ) + fpj+1 − fpj , (33)

4.4.9

Analysis of Stationary Discrete Shocks

where u ¯ and w ¯ are the arithmetic averages u ¯=

The ideal model of a discrete shock is illustrated in ﬁgure (4). Suppose that wL and wR are left and right states which satisfy the jump conditions for a stationary shock, and that the corresponding ﬂuxes are fL = f (wL ) and fR = f (wR ). Since the shock is stationary fL = fR . The ideal discrete shock has constant states wL to the left and wR to the right, and a single point with an intermediate value wA . The intermediate value is needed to allow the discrete solution to correspond to a true solution in which the shock wave does not coincide with an interface between two mesh cells.

1 1 ¯ = (wj+1 + wj ) . (uj+1 + uj ) , w 2 2

Thus these schemes are closely related to schemes which introduce separate splittings of the convective and pressure terms, such as the wave-particle scheme [141, 8], the advection upwind splitting method (AUSM) [106, 175], 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 ﬁrst order diﬀusive ﬂux of the form dj+ 12

1 = αj+ 12 Bj+ 12 (wj+1 − wj ) , 2

Schemes corresponding to one, two or three terms in equation (35) are examined in [80]. The analysis of these three cases shows that a discrete shock structure with a single interior point is supported by artiﬁcial diﬀusion that satisﬁes the two conditions that

(34)

where the matrix Bj+ 12 determines the properties of the scheme and the scaling factor αj+ 12 is included for convenience. All the previous schemes can be obtained by representing Bj+ 12 as a polynomial in the matrix Aj+ 12 deﬁned by equation (26). Schemes of this class were considered by Van Leer [99]. According to the Cayley-Hamilton theorem, a matrix satisﬁes its own characteristic equation. Therefore the third and higher powers of A can be eliminated, and there is no loss of generality in limiting Bj+ 12 to a polynomial of degree 2, Bj+ 12 = α0 I + α1 Aj+ 12 + α2 A2j+ 1 . 2

1. it produces an upwind ﬂux if the ﬂow is determined to be supersonic through the interface 2. it satisﬁes a generalized eigenvalue problem for the exit from the shock of the form (AAR − αAR BAR ) (wR − wA ) = 0,

(36)

where AAR is the linearized Jacobian matrix and BAR is the matrix deﬁning the diﬀusion for the interface AR. This follows from the equilibrium condition hRA = hRR for the cell j +1 in ﬁgure 4. These two conditions are satisﬁed by both the characteristic scheme and also the CUSP scheme, provided that the coeﬃcients of convective diﬀusion and pressure diﬀerences are correctly balanced. Scalar diﬀusion does not satisfy the ﬁrst condition. In the case of the CUSP scheme (30) equation 36 reduces to α∗ c (wR − wA ) = 0 ARA + 1+β

(35)

The characteristic upwind scheme for which Bj+ 12 = Aj+ 12 is obtained by substituting Aj+ 12 = T ΛT −1,

A2j+ 1 = T Λ2 T −1 . Then α0 , α1 , and α2 are deter2 mined from the three equations α0 + α1 λk + α2 λ2k = |λk | , k = 1, 2, 3.

The same representation remains valid for three dimensional ﬂow because Aj+ 12 still has only three distinct eigenvalues u, u + c, u − c. 12

Thus wR − wA is an eigenvector of the Roe maα∗ c trix ARA , and − 1+β is the corresponding eigenvalue. Since the eigenvalues are u, u + c, and u − c, the only choice which leads to positive diﬀusion when u > 0 is u − c, yielding the relationship

the discrepancy between the convective terms ρ u ρu , ρH in the ﬂux vector, which contain ρH, and the state vector which contains ρE. This may be remedied by introducing a modiﬁed state vector ρ wh = ρu . ρH

α∗ c = (1 + β)(c − u), 0 < u < c Thus there is a one parameter family of schemes which support the ideal shock structure. The term β(fR − fA ) contributes to the diﬀusion of the convective terms. Allowing for the split (31), the total eﬀective coeﬃcient of convective diﬀusion is αc = α∗ c + βu. A CUSP scheme with low numerical diffusion is then obtained by taking α = |M |, leading to the coeﬃcients illustrated in ﬁgure 5

α(M)

Then one introduces the linearization fR − fL = Ah (whR − whL ). Here Ah may be calculated in the same way as the standard Roe linearization. Introduce the weighted averages deﬁned by equation (27). Then

β (M)

1

1

0

u Ah = − γ+1 γ 2 −uH

2

-1

1 M

-1

1

M

-1

0

γ+1 γ u

γ−1 γ

H

u

.

The eigenvalues of Ah are u, λ+ and λ− where γ + 1 γ + 1 2 c2 − u 2 u± ( u) + . (37) λ± = 2γ 2γ γ

Figure 5: Diﬀusion Coeﬃcients.

4.4.10

1

Now both CUSP and characteristic schemes which preserve constant stagnation enthalpy in steady ﬂow can be constructed from the modiﬁed Jacobian matrix Ah [80]. These schemes also produce a discrete shock structure with one interior point in steady ﬂow. Then one arrives at four variations with this property, which can conveniently be distinguished as the E- and H-CUSP schemes, and the E- and Hcharacteristic schemes.

CUSP and Characteristic Schemes Admitting Constant Total Enthalpy in Steady Flow

In steady ﬂow 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 semi-discrete schemes do not necessarily satisfy this property. In the case of a semidiscrete scheme expressed in viscosity form, equations (11) 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 ρ replaced by ρH. When the standard characteristic decomposition (28) is used, the viscous ﬂuxes for ρ and ρH which result from composition of the ﬂuxes for the characteristic variables do not have this property, and H is not constant in the discrete solution. In practice there is an excursion of H in the discrete shock structure which represents a local heat source. In very high speed ﬂows the corresponding error in the temperature may lead to a wrong prediction of associated eﬀects such as chemical reactions.

4.5

Multidimensional Schemes

The simplest approach to the treatment of multidimensional problems on structured meshes is to apply the one-dimensional 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 ﬂow rather than the mesh. A thorough review is given by Pailliere and DeConinck in reference [132]. Residual distribution schemes are an attractive approach for triangulated meshes. In these the residual deﬁned by the space derivatives is evaluated for

The source of the error in the stagnation enthalpy is 13

each cell, and then distributed 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 diﬀusion in the direction normal to the ﬂow. For the Euler equations the residual can be linearized by assuming that the parameter vector with com√ √ √ ponents ρ, ρui , and ρH varies linearly over the cell. Then ∂w ∂fj (w) = Aj ∂xj ∂xj

[42, 36, 45, 136, 180].

4.6

of

the

Viscous

The discretization of the viscous terms of the Navier Stokes equations requires an approximation to the ∂ui in order to calculate the tenvelocity derivatives ∂x j sor σij , equation (3). Then the viscous terms may be included in the ﬂux balance (4). In order to evaluate the derivatives one may apply the Gauss formula to a control volume V with the boundary S. ∂ui dv = ui nj ds V ∂xj S

∂f

where the Jacobian matrices Aj = ∂wj are evaluated with Roe averaging of the values of w at the vertices. Waves in the direction n can then be expressed in terms of the eigenvectors of nj Aj , and a positive distribution scheme is used for waves in preferred directions. The best choice of these directions is the subject of ongoing research, but preliminary results indicate the possibility of achieving high resolution of shocks and contact discontinuities which are not aligned with mesh lines [132].

where nj is the outward normal. For a tetrahedral or hexahedral cell this gives ∂ui 1 = u i nj s ∂xj vol faces

Hirsch and Van Ransbeeck adopt an alternative approach in which they directly construct directional diﬀusive terms on structured meshes, with antidiﬀusion controlled by limiters based on comparisons of slopes in diﬀerent directions [60]. They also show promising results in calculations of nozzles with multiply reﬂected oblique shocks.

4.5.1

Discretization Terms

(38)

where ui is an estimate of the average of ui over the face. If u varies linearly over a tetrahedral cell this is exact. Alternatively, assuming a local transformation to computational coordinates ξj , one may apply the chain rule −1 ∂u ∂u ∂ξ ∂u ∂x = (39) = ∂x ∂ξ ∂x ∂ξ ∂ξ

High Order Godunov Schemes, and Kinetic Flux Splitting

i Here the transformation derivatives ∂x ∂ξj can be evaluated by the same ﬁnite diﬀerence formulas as the ∂u i velocity derivatives ∂u ∂ξj In this case ∂ξ is exact if u is a linearly varying function.

A substantial body of current research is directed toward the implementation of truly multi-dimensional upwind schemes [59, 135, 101]. Reference [132] provides a thorough review of recent developments in this ﬁeld. Some of the most impressive simulations of time dependent ﬂows with strong shock waves have been achieved with higher order Godunov schemes [179]. In these schemes the average value in each cell is updated by applying the integral conservation law using interface ﬂuxes predicted from the exact or approximate solution of a Riemann problem between adjacent cells. A higher order estimate of the solution is then reconstructed from the cell averages, and slope limiters are applied to the reconstruction. An example is the class of essentially non-oscillatory (ENO) schemes, which can attain a very high order of accuracy at the cost of a substantial increase in computational complexity [32, 153, 151, 152]. Methods based on reconstruction can also be implemented on unstructured meshes [13, 12]. Recently there has been an increasing interest in kinetic ﬂux splitting schemes, which use solutions of the Boltzmann equation or the BGK equation to predict the interface ﬂuxes

For a cell-centered discretization (ﬁgure 6a) ∂u ∂ξ is needed at each face. The simplest procedure is to ∂u evaluate ∂u ∂ξ in each cell, and to average ∂ξ between the two cells on either side of a face [87]. The resulting discretization does not have a compact stencil, and supports undamped oscillatory modes. In a one 2 dimensional calculation, for example, ∂∂xu2 would be i +ui−2 . In order to produce a discretized as ui+2 −2u 4∆x2 compact stencil ∂u may be estimated from a control ∂x volume centered on each face, using formulas (38) or (39) [144]. This is computationally expensive because the number of faces is much larger than the number of cells. In a hexahedral mesh with a large number of vertices the number of faces approaches three times the number of cells. This motivates the introduction of dual meshes for the evaluation of the velocity derivatives and the ﬂux balance as sketched in ﬁgure 6. The ﬁgure shows both cell-centered and cell-vertex schemes. The dual mesh connects cell centers of the primary mesh. If there is a kink in the primary mesh, the dual cells 14

dual cell

stituted for the solution at the vertices appearing in the local stencil. Figure 7 illustrates the discretization of the Laplacian uxx + uyy which is obtained with linear elements. It shows a particular triangulation such that the approximation is locally consistent with uxx + 3uyy . Thus the use of an irregular triangulation in the boundary layer may signiﬁcantly degrade the accuracy.

σi j

3

6a: Cell-centered scheme. 6b: Cell-vertex scheme. σij evaluated at vertices σij evaluated at cell cenof the primary mesh ters of the primary mesh

d

h

-6

Figure 6: Viscous discretizations for cell-centered and cell-vertex algorithms.

Coefficients resulting from linear elements

e

h

1 a

should be formed by assembling contiguous fractions of the neighboring primary cells. On smooth meshes comparable results are obtained by either of these formulations [114, 115, 107]. If the mesh has a kink the cell-vertex scheme has the advantage that the ∂ui derivatives ∂x are calculated in the interior of a j regular cell, with no loss of accuracy.

1c

1 b h

h

Figure 7: Example of discretization uxx + uyy on a triangular mesh. The discretization is locally equivb +uc , 3uyy = alent to the approximation uxx = ua −2u h2 3ud −6ue +3ub . h2

A desirable property is that a linearly varying velocity distribution, as in a Couette ﬂow, should produce a constant stress and hence an exact stress balance. This property is not necessarily satisﬁed in general by ﬁnite diﬀerence or ﬁnite volume schemes on curvilinear meshes. The characterization k-exact has been proposed for schemes that are exact for polynomials of degree k. The cell-vertex ﬁnite volume scheme is linearly exact if the derivatives are ∂ui is exevaluated by equation (39), since then ∂x j actly evaluated as a constant, leading to constant viscous stresses σij , 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 ﬁgure 6 still yields a perfect balance. The use of ∂ui , however, requires the equation (39) to evaluate ∂x j additional calculation or storage of the nine metric ∂ui quantities ∂x in each cell, whereas equation (38) j can be evaluated from the same face areas that are used for the ﬂux balance.

4.7

Time Stepping Schemes

If the space discretization procedure is implemented separately, it leads to a set of coupled ordinary differential equations, which can be written in the form dw + R(w) = 0, dt

(40)

where w is the vector of the ﬂow variables at the mesh points, and R(w) is the vector of the residuals, consisting of the ﬂux balances deﬁned by the space discretization scheme, together with the added dissipative terms. If the objective is simply to reach the steady state and details of the transient solution are immaterial, the time-stepping scheme may be designed solely to maximize the rate of convergence. The ﬁrst decision that must be made is whether to use an explicit scheme, in which the space derivatives are calculated from known values of the ﬂow variables at the beginning of the time step, or an implicit scheme, in which the formulas for the space derivatives include as yet unknown values of the ﬂow variables at the end of the time step, leading to the need to solve coupled equations for the new values. The permissible time step for an explicit scheme is limited by the Courant-Friedrichs-Lewy (CFL) condition, which states that a diﬀerence scheme cannot be a convergent and stable approximation unless its

In the case of an unstructured mesh, the weak form (6) leads to a natural discretization with linear elements, in which the piecewise linear approximation yields a constant stress in each cell. This method yields a representation which is globally correct when averaged over the cells, as is proved by energy estimates for elliptic problems [15]. It should be noted, however, that it yields formulas that are not necessarily locally consistent with the diﬀerential equations, if Taylor series expansions are sub15

domain of dependence contains the domain of dependence of the corresponding diﬀerential equation. One can anticipate that implicit schemes will yield convergence in a smaller number of time steps, because the time step is no longer constrained by the CFL condition. Implicit schemes will be eﬃcient, however, only if the decrease in the number of time steps outweighs the increase in the computational effort per time step consequent upon the need to solve coupled equations. The prototype implicit scheme w at t + µ∆t as can be formulated by estimating ∂∂t n a linear combination of R(w ) and R(wn+1 ). The resulting equation wn+1 = wn − ∆t (1 − µ) R (wn ) + µR wn+1

scheme should have a stability region which contains a substantial interval of the negative real axis, as well as an interval along the imaginary axis. To achieve this it pays to treat the convective and dissipative terms in a distinct fashion. Thus the residual is split as R(w) = Q(w) + D(w), where Q(w) is the convective part and D(w) the dissipative part. Denote the time level n∆t by a superscript n. Then the multistage time stepping scheme is formulated as w(n+1,0) w(n+1,k)

can be linearized as ∂R I + µ∆t δw + ∆tR (wn ) = 0. ∂w

wn+1

= wn ...

= wn − αk ∆t Q(k−1) + D(k−1) ... =

w(n+1,m) ,

where the superscript k denotes the k-th stage, αm = 1, and

If one sets µ = 1 and lets ∆t → ∞ this reduces to the Newton iteration , which has been successfully used in two-dimensional calculations [172, 50]. In the three-dimensional case with, say, an N × N × N mesh, the bandwidth of the matrix that must be inverted is of order N 2 . Direct inversion requires a number of operations proportional to the number of unknowns multiplied by the square of the bandwidth of the order of N 7 . This is prohibitive, and forces recourse to either an approximate factorization method or an iterative solution method.

Q(0) Q(k) D(k)

=

Q (wn ) , D(0) = D (wn )

...

= Q w(n+1,k)

= βk D w(n+1,k) + (1 − βk )D(k−1) .

The coeﬃcients αk are chosen to maximize the stability interval along the imaginary axis, and the coeﬃcients βk are chosen to increase the stability interval along the negative real axis.

Alternating direction methods, which introduce factors corresponding to each coordinate, are widely used for structured meshes [17, 137]. They cannot be implemented on unstructured tetrahedral meshes that do not contain identiﬁable mesh directions, although other decompositions are possible [108]. If one chooses to adopt the iterative solution technique, the principal alternatives are variants of the Gauss-Seidel and Jacobi methods. A symmetric Gauss-Seidel method with one iteration per time step is essentially equivalent to an approximate lower-upper (LU) factorization of the implicit scheme [86, 125, 31, 183]. On the other hand, the Jacobi method with a ﬁxed number of iterations per time step reduces to a multistage explicit scheme, belonging to the general class of RungeKutta schemes [33]. Schemes of this type have proved very eﬀective for wide variety of problems, and they have the advantage that they can be applied equally easily on both structured and unstructured meshes [84, 67, 69, 145].

These schemes do not fall within the standard framework of Runge-Kutta schemes, and they have much larger stability regions [69]. Two schemes which have been found to be particularly eﬀective are tabulated below. The ﬁrst is a four-stage scheme with two evaluations of dissipation. Its coeﬃcients are α1 α2 α3 α4

= 13 4 = 15 5 = 9 =1

β1 β2 β3 β4

=1 = 12 . =0 =0

(41)

The second is a ﬁve-stage scheme with three evaluations of dissipation. Its coeﬃcients are α1 α2 α3 α4 α5

4.8

If one reduces the linear model problem corresponding to (40) to an ordinary diﬀerential equation by substituting a Fourier mode w ˆ = eipxj , the resulting Fourier symbol has an imaginary part proportional to the wave speed, and a negative real part proportional to the diﬀusion. Thus the time stepping

4.8.1

= 14 = 16 = 38 = 12 =1

β1 β2 β3 β4 β5

=1 =0 = 0.56 . =0 = 0.44

(42)

Multigrid Methods Acceleration of Steady Flow Calculations

Radical improvements in the rate of convergence to a steady state can be realized by the multigrid time16

stepping technique. The concept of acceleration by the introduction of multiple grids was ﬁrst proposed by Fedorenko [48]. There is by now a fairly welldeveloped theory of multigrid methods for elliptic equations based on the concept that the updating scheme acting as a smoothing operator on each grid [24, 53]. This theory does not hold for hyperbolic systems. Nevertheless, it seems that it ought to be possible to accelerate the evolution of a hyperbolic system to a steady state by using large time steps on coarse grids so that disturbances will be more rapidly expelled through the outer boundary. Various multigrid time-stepping schemes designed to take advantage of this eﬀect have been proposed [124, 65, 55, 71, 29, 6, 57, 83, 93].

E

T

E

8a: 3 Levels. E

E

E

E

E

E

T

E

T

E

E

E

E

T

E

8b: 4 Levels. E

E

E

E

4 Level Cycle

T

4 Level Cycle

8c: 5 Levels. Figure 8: Multigrid W -cycle for managing the grid calculation. E, evaluate the change in the ﬂow for one step; T , transfer the data without updating the solution.

(0)

wk = Tk,k−1 wk−1 , where wk−1 is the current value on grid k − 1, and Tk,k−1 is a transfer operator. Next it is necessary to transfer a residual forcing function such that the solution grid k is driven by the residuals calculated on grid k − 1. This can be accomplished by setting (0) Pk = Qk,k−1 Rk−1 (wk−1 ) − Rk wk ,

and consequently the very large eﬀective time step of the complete cycle costs only slightly more than a single time step in the ﬁne grid.

where Qk,k−1 is another transfer operator. Then Rk (wk ) is replaced by Rk (wk )+Pk in the time- stepping scheme. Thus, the multistage scheme is reformulated as (1) (0) (0) = wk − α1 ∆tk Rk + Pk wk ···

E

E

One can devise a multigrid scheme using a sequence of independently generated coarser meshes by eliminating alternate points in each coordinate direction. In order to give a precise description of the multigrid scheme, subscripts may be used to indicate the grid. Several transfer operations need to be deﬁned. First the solution vector on grid k must be initialized as

(q+1) wk

E

E

4.8.2

Multigrid Implicit Schemes for Unsteady Flow

Time dependent calculations are needed for a number of important applications, such as ﬂutter analysis, or the analysis of the ﬂow past a helicopter rotor, in which the stability limit of an 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].

··· (0) (q) = wk − αq+1 ∆tk Rk + Pk . (m)

The result wk then provides the initial data for grid k + 1. Finally, the accumulated correction on grid k has to be transferred back to grid k − 1 with the aid of an interpolation operator Ik−1,k . With properly optimized coeﬃcients multistage time-stepping schemes can be very eﬃcient drivers of the multigrid process. A W -cycle of the type illustrated in Figure 8 proves to be a particularly eﬀective strategy for managing the work split between the meshes. In a three-dimensional case the number of cells is reduced by a factor of eight on each coarser grid. On examination of the ﬁgure, it can therefore be seen that the work measured in units corresponding to a step on the ﬁne grid is of the order of

Suppose that (40) is approximated as Dt wn+1 + R(wn+1 ) = 0. Here Dt is a k th order accurate backward diﬀerence operator of the form 1 1 − q (∆ ) , ∆t q=1 q k

Dt = where

∆− wn+1 = wn+1 − wn .

1 + 2/8 + 4/64 + . . . < 4/3, 17

Applied to the linear diﬀerential equation

nates this has the form τ 2 − βτ M β2 M τ − M τ2 + 1 β β P = 0 0 0 0 0 0

dw = αw dt the schemes with k = 1, 2 are stable for all α∆t in the left half plane (A-stable). Dahlquist has shown that A-stable linear multi-step schemes are at best second order accurate [38]. Gear however, has shown that the schemes with k ≤ 6 are stiﬄy stable [49], and one of the higher order schemes may oﬀer a better compromise between accuracy and stability, depending on the application.

0 0 0 0 1

$ = τ = 1 − M 2 , if M < 1 % $ 1 1 − M 2 , τ = 1 − 2 , if M ≥ 1 β = M Turkel has proposed an asymmetric preconditioner which has also proved eﬀective, particularly for ﬂow at low Mach numbers [171]. The use of these preconditioners can lead to instability at stagnation points where there is a zero eigenvalue which cannot be equalized with the eigenvalues ±c. β

The preconditioners of Van Leer and Turkel do not take account of the eﬀect of diﬀerences in the mesh intervals in the diﬀerent coordinate directions. The need to resolve the boundary layer generally compels the introduction of mesh cells with very high aspect ratios near the boundary, and these can lead to a severe reduction in the rate of convergence to a steady state. Pierce has recently obtained impressive results using diagonal and block-Jacobi preconditioners which include the mesh intervals [133].

∂w = R∗ (w), ∂t∗ where 2 n 3 1 n−1 w + R(w) + w − w , 2∆t ∆t 2∆t

and the last two terms are treated as ﬁxed source terms. The ﬁrst term shifts the Fourier symbol of the equivalent model problem to the left in the complex plane. While this promotes stability, it may also require a limit to be imposed on the magnitude of the local time step ∆t∗ relative to that of the implicit time step ∆t. This may be relieved by a pointimplicit modiﬁcation of the multi-stage scheme [119]. In the case of problems with moving boundaries the equations must be modiﬁed to allow for movement and deformation of the mesh.

An alternative approach has recently been proposed by Ta’asan [167], in which the equations are written in a canonical form which separates the equations describing acoustic waves from those describing convection. In terms of the velocity components u, v and the vorticity ω, temperature T , entropy s and total enthalpy H, the equations describing steady two-dimensional ﬂow can be written as 0 − D1 D2 0 u ∂ ∂ − ∂y ∂x −1 0 0 v c2 1 0 0 −q − γ(γ−1) D3 q D3 ω = 0 0 s 0 0 T ρQ 0 H 0 0 0 0 ρQ

This method has proved eﬀective for the calculation of unsteady ﬂows that might be associated with wing ﬂutter [3, 4] and also in the calculation of unsteady incompressible ﬂows [18]. It has the advantage that it can be added as an option to a computer program which uses an explicit multigrid scheme, allowing it to be used for the eﬃcient calculation of both steady and unsteady ﬂows.

4.9

0 0 0 τ 0

where

Equation (40) is now treated as a modiﬁed steady state problem to be solved by a multigrid scheme using variable local time steps in a ﬁctitious time t∗ . For example, in the case k = 2 one solves

R∗ (w) =

0 0 τ 0 0

where D1

=

D2

=

D3

=

Preconditioning

Another way to improve the rate of convergence to a steady state is to multiply the space derivatives in equation (1) by a preconditioning matrix P which is designed to equalize the eigenvalues, so that all the waves can be advanced with optimal time steps. A symmetric preconditioner which equalizes the eigenvalues has been proposed by Van Leer [102]. When the equations are written in stream-aligned coordi-

Q =

∂ ρ 2 2 ∂ − uv (c − u ) c2 ∂x ∂y ∂ ρ 2 2 ∂ − uv (c − u ) c2 ∂y ∂x ∂ ∂ −u v ∂x ∂y ∂ ∂ u +v ∂x ∂y

Here the ﬁrst two equations describe an elliptic system if the ﬂow is subsonic, while the remaining equations are convective. Now separately optimized multigrid procedures are used to solve the two sets of equations, which are essentially decoupled. 18

4.10

5

High Order Schemes and Mesh Reﬁnement

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 improving the accuracy are to increase the order of the discrete scheme, and reduce the mesh interval. High order diﬀerence methods are most easily implemented on Cartesian, or at least extremely smooth grids. The expansion of the stencil as the order is increased leads to the need for complex boundary conditions. Compact schemes keep the stencil as small as possible [140, 104, 28]. On simple domains, spectral methods are particularly eﬀective, especially in the case of periodic boundary conditions, and can be used to produce exponentially fast convergence of the error as the mesh interval is decreased [127, 27]. A compromise is to divide the ﬁeld into subdomains and introduce high order elements. This approach is used in the spectral element method [92].

CURRENT STATUS OF NUMERICAL SIMULATION

This section presents some representative numerical results which conﬁrm the properties of the algorithms which have been reviewed in the last section. These have been drawn from the work of the author and his associates. They also illustrate the kind of calculation which can be performed in an industrial environment, where rapid turn around is important to allow the quick assessment of design changes, and computational costs must be limited.

5.1

One dimensional shock

In order to verify the discrete structure of stationary shocks, calculations were performed for a one dimensional problem with initial data containing left and right states compatible with the Rankine Hugoniot conditions. An intermediate state consisting of the arithmetic average of the left and right states was introduced at a single cell in the center of the domain. With this intermediate state the system is not in equilibrium, and the time dependent equations were solved to ﬁnd an equilibrium solution with a stationary shock wave separating the left and right states. Table 1 shows the result for a shock wave at Mach 20. This calculation used the H-CUSP scheme, which allows a solution with constant stagnation enthalpy. The SLIP-JST construction was used with the limiter deﬁned by equation (23), and q = 3. The table shows the values of ρ, u, H, p, M and the entropy p − log pρLγ . A perfect one point shock S = log ργ L structure is displayed. The entropy is zero to 4 decimal places upstream of the shock, exhibits a slight excursion at the interior point, and is constant to 4 decimal places downstream of the shock. It may be noted that the mass, momentum and energy of the initial data are not compatible with the ﬁnal equilibrium state. According to conservation arguements the total mass, momentum and energy must remain constant if the outﬂow ﬂux fR remains equal to the inﬂow ﬂux fL . Therefore fR must be allowed to vary according to an appropriate outﬂow boundary condition to allow the total mass, momentum and energy to be adjusted to values compatible with equilibrium.

High order diﬀerence schemes and spectral methods have proven particularly useful in direct NavierStokes simulations of transient and turbulent ﬂows. High order methods are also beneﬁcial in computational aero-acoustics, where it is desired to track waves over long distances with minimum error. If the ﬂow contains shock waves or contact discontinuities, the ENO method may be used to construct high order non-oscillatory schemes. In multi-dimensional ﬂow simulations, global reduction of the mesh interval can be prohibitively expensive, motivating the use of adaptive mesh reﬁnement procedures which reduce the local mesh width h if there is an indication that the error is too large [21, 39, 109, 61, 138, 103]. In such h-reﬁnement methods, simple error indicators such as local solution gradients may be used. Alternatively, the discretization error may be estimated by comparing quantities calculated with two mesh widths, 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 calculation is completed. This kind of local adaptive control can also be applied to the local order of a ﬁnite element method to produce a p-reﬁnement method, where p represents the order of the polynomial basis functions. Finally, both h- and p- reﬁnement can be combined to produce an h-p method in which h and p are locally optimized to yield a solution with minimum error [126]. Such methods can achieve exponentially fast convergence, and are well established in computational solid mechanics.

5.2

Euler Calculations for Airfoils and Wings

The results of transonic ﬂow calculations for two well known airfoils, the RAE 2822 and the NACA 19

I 19 20 21 22 23 24 25

ρ 1.0000 1.0000 1.0000 4.1924 5.9259 5.9259 5.9259

H 283.5000 283.5000 283.5000 283.4960 283.4960 283.4960 283.4960

p 1.0000 1.0000 1.0000 307.4467 466.4889 466.4889 466.4889

M 20.0000 20.0000 20.0000 0.7229 0.3804 0.3804 0.3804

a 160x32 mesh. Mesh

40x8 80x16 160x32

Table 1: Shock Wave at Mach 20

RAE 2822 Mach .50 α 3◦ .0062 .0013 .0000

NACA 0012 Mach .50 α 3◦ .0047 .0008 .0000

Korn Airfoil Mach .75 α 0◦ .0098 .0017 .0000

Table 2: Drag Coeﬃcient on a sequence of meshes

0012, are presented in ﬁgures (19-21). The H-CUSP scheme was again used with the SLIP-JST construction. The limiter deﬁned by equation (23) was used with q = 3. The 5 stage time stepping scheme (42) was augmented by the multigrid scheme described in section 4.2 to accelerate convergence to a steady state. The equations were discretized on meshes with O-topology extending out to a radius of about 100 chords. In each case the calculations were performed on a sequence of successively ﬁner meshes from 40x8 to 320x64 cells, while the multigrid cycles on each of these meshes descended to a coarsest mesh of 10x2 cells. Figures 19-21 show the ﬁnal results on 320x64 meshes for the RAE 2822 airfoil at Mach .75 and 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 coeﬃcient Cp = 1p−p 2 2 ρ∞ q∞ is plotted with the negative (suction) pressures upward, so that the upper curve represents the ﬂow over the upper side of a lifting airfoil. The convergence histories show the mean rate of change of the density, and also the total number of supersonic points in the ﬂow ﬁeld, which provides a useful measure of the global convergence of transonic ﬂow calculations such as these. In each case the convergence history is shown for 100 cycles, while the pressure distribution is displayed after a suﬃcient number of cycles for its convergence. The pressure distribution of the RAE 2822 airfoil converged in only 25 cycles. Convergence was slower for the NACA 0012 airfoil. In the case of ﬂow at Mach .8 and 1.25◦ angle of attack, additional cycles were needed to damp out a wave downstream of the weak shock wave on the lower surface.

As a further test of the performance of the H-CUSP scheme, the ﬂow past the ONERA M6 wing was calculated on a mesh with C-H topology and 192x32x48 = 294912 cells. Figure 22 shows the result at Mach .84 and 3.06◦ angle of attack. This again veriﬁes the non-oscillatory character of the solution, and the sharp resolution of shock waves. In this case 50 cycles were suﬃcient for convergence of the pressure distributions. Figure 9 shows a calculation of the Northrop YF23 by R.J. Busch, Jr., who used the author’s FLO57 code to solve the Euler equations [26]. Although an inviscid model of the ﬂow was used, it can be seen that the simulations are in good agreement with wind tunnel measurements both at Mach .90, with angles of attack of 0, 8 and 16 degrees, and at Mach 1.5 with angles of attack of 0, 4 and 8 degrees. At a high angle of attack the ﬂow separates from the leading edge, and this example shows that in situations where the point of separation is ﬁxed, an inviscid model may still produce a useful prediction. Thus valuable information for the aerodynamic design could be obtained with a relatively inexpensive computational model. The next ﬁgures show the results of calculations using the AIRPLANE code developed by T.J. Baker and the author, to solve the Euler equations on an unstructured mesh. This provides the ﬂexibility to treat arbitrarily complex conﬁgurations without the need to spend months developing an acceptable mesh. Figures 10 and 11 show calculations for supersonic transport conﬁgurations which were performed by Susan Cliﬀ. The agreement with experimental data is quite good, and it has also been possible to predict the sonic boom signature [34]. Figure 12 shows an Euler calculation for the McDonnell Douglas MD11 with ﬂow through the engine nacelles, using 348407 mesh points of 2100466 tetrahedra. This calculation takes 4 hours on an IBM 590 workstation. A parallel version of the code has been developed in collaboration with W.S. Chen, and the same calculation can be performed in 16 minutes using 16 processors of an IBM SP2. The parallel speed-up for the MD11 is shown in table 3.

As a further check on accuracy the drag coeﬃcient should be zero in subsonic ﬂow, or in shock free transonic ﬂow. Table 2 shows the computed drag coeﬃcient on a sequence of three meshes for three examples. The ﬁrst two are subsonic ﬂows over the RAE 2822 and NACA 0012 airfoils at Mach .5 and 3◦ angle of attack. The third is the ﬂow over the shock free Korn airfoil at its design point of Mach .75 and 0◦ angle of attack. In all three cases the drag coeﬃcient is calculated to be zero to four digits on 20

Figure 10: Comparison of Experimental and Calculated Results for a HSCT Conﬁguration

Figure 9: Comparison of Experimental and Computed Drag Rise Curve for the YF-23 (Supplied by R. J. Bush Jr.) No. of Nodes 1 2 4 8 16

Seconds/Cycle 36.03 18.11 9.11 4.66 2.39

Speedup 1.00 1.99 3.96 7.73 15.08

both the experimental data and the visualization are quite good. Figure 14 shows an unsteady ﬂow calculation for a pitching airfoil performed by J. Alonso using the code UFLO82, which he jointly developed with L. Martinelli and the author [4]. This uses the multigrid implicit scheme described in Section 3.7.2 which allows the number of time steps to be reduced from several thousand to 36 per pitching cycle. The agreement with experimental data is quite good.

Table 3: AIRPLANE Parallel Performance on the SP2, MD-11 Model

5.3

Viscous Flow Calculations

The next ﬁgures show viscous simulations based on the solution of the Reynolds averaged Navier Stokes equations with turbulence models. Figure 13 shows a two-dimensional calculation for the RAE 2822 airfoil by L. Martinelli. The vertical axis represents the negative pressure coeﬃcient, and there is a shock wave half way along the upper surface. This example conﬁrms that in the absence of signiﬁcant shock induced separation, simulations performed on a sufﬁciently ﬁne mesh (512 x 64) can produce excellent agreement with experimental data. Figure 18 shows a simulation of the McDonnell-Douglas F18 performed by R.M. Cummings, Y.M. Rizk, L.B. Schiﬀ and N.M. Chaderjian at NASA Ames [37]. They used a multiblock mesh with about 900000 mesh points. While this is probably not enough for an accurate quantitative prediction, the agreement with

5.4

Ship Wave Resistance calculations

Figures 51-53 show the results of an application of the same multigrid ﬁnite volume techniques to the calculation of the ﬂow past a naval frigate, using a code which was developed by J. Farmer, L. Martinelli 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 non-linear boundary condition, while artiﬁcial compressibility was used to reat the incompressible ﬂow equations. 21

Figure 11: Pressure Contours and Sonic Boom on a Representative HSCT Conﬁguration

Figure 13: Two-Dimensional Turbulent Viscous Calculation (by Luigi Martinelli) improvements in aerodynamic eﬃciency. The simplest approach to optimization is to deﬁne the geometry through a set of design parameters, which may, for example, be the weights αi applied to a set of shape functions bi (x) so that the shape is represented as

Figure 12: Computed Pressure Field for a McDonnell Douglas MD11

6

6.1

AERODYNAMIC OPTIMIZATION

SHAPE

f (x) =

αi bi (x).

Then a cost function I is selected which might, for example, be the drag coeﬃcient or the lift to drag ratio, and I is regarded as a function of the parame∂I may now be estimated ters αi . The sensitivities ∂α i by making a small variation δαi in each design parameter in turn and recalculating the ﬂow to obtain the change in I. Then

Optimization and Design

Traditionally the process of selecting design variations has been carried out by trial and error, relying on the intuition and experience of the designer. With currently available equipment the turn around for numerical simulations is becoming so rapid that it is feasible to examine an extremely large number of variations. It is not at all likely that repeated trials in an interactive design and analysis procedure can lead to a truly optimum design. In order to take full advantage of the possibility of examining a large design space the numerical simulations need to be combined with automatic search and optimization procedures. This can lead to automatic design methods which will fully realize the potential

I(αi + δαi ) − I(αi ) ∂I ≈ . ∂αi δαi ∂I The gradient vector ∂α may now be used to determine a direction of improvement. The simplest procedure is to make a step in the negative gradient direction by setting

αn+1 = αn − λδα, 22

inverse problem as a special case of the optimization problem, with a cost function which measures the error in 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 error, 1 I= (p − pd )2 dB, 2 B or possibly a more general Sobolev norm of the pressure error. This has the advantage of converting a possibly ill posed problem into a well posed one. It has the disadvantage that it incurs the computational costs associated with optimization procedures.

6.2

In order to reduce the computational costs, it turns out that there are advantages in formulating both the inverse problem and more general aerodynamic problems within the framework of the mathematical theory for the control of systems governed by partial diﬀerential equations [105]. A wing, for example, is a device to produce lift by controlling the ﬂow, and its design can be regarded as a problem in the optimal control of the ﬂow equations by variation of the shape of the boundary. If the boundary shape is regarded as arbitrary within some requirements of smoothness, then the full generality of shapes cannot be deﬁned with a ﬁnite number of parameters, and one must use the concept of the Frechet derivative of the cost with respect to a function. Clearly, such a derivative cannot be determined directly by ﬁnite diﬀerences of the design parameters because there are now an inﬁnite number of these. Using techniques of control theory, however, the gradient can be determined indirectly by solving an adjoint equation which has coeﬃcients deﬁned by the solution of the ﬂow equations. The cost of solving the adjoint equation is comparable to that of solving the ﬂow equations. Thus the gradient can be determined with roughly the computational costs of two ﬂow solutions, independently of the number of design variables, which may be inﬁnite if the boundary is regarded as a free surface.

Figure 14: Mach Number Contours. Pitching Airfoil Case. Re = 1.0 × 106 , M∞ = 0.796, Kc = 0.202. so that to ﬁrst order I + δI = I −

Application of Control Theory

∂I T ∂I ∂I T δα = I − λ . ∂α ∂α ∂α

More sophisticated search procedures may be used such as quasi-Newton methods, which attempt to es2 I timate the second derivative ∂α∂i ∂α of the cost funcj ∂I in successive tion from changes in the gradient ∂α optimization steps. These methods also generally introduce line searches to ﬁnd the minimum in the search direction which is deﬁned at each step. The main disadvantage of this approach is the need for a number of ﬂow calculations proportional to the number of design variables to estimate the gradient. The computational costs can thus become prohibitive as the number of design variables is increased.

An alternative approach is to cast the design problem as a search for the shape that will generate the desired pressure distribution. This approach recognizes that the designer usually has an idea of the the kind of pressure distribution that will lead to the desired performance. Thus, it is useful to consider the inverse problem of calculating the shape that will lead to a given pressure distribution. The method has the advantage that only one ﬂow solution is required to obtain the desired design. Unfortunately, a physically realizable shape may not necessarily exist, unless the pressure distribution satisﬁes certain constraints. Thus the problem must be very carefully formulated; otherwise it may be ill posed.

For ﬂow about an airfoil or wing, the aerodynamic properties which deﬁne the cost function are functions of the ﬂow-ﬁeld variables (w) and the physical location of the boundary, which may be represented by the function F , say. Then I = I (w, F ) , and a change in F results in a change

The diﬃculty that the target pressure may be unattainable may be circumvented by treating the

δI = 23

∂I T ∂I T δw + δF , ∂w ∂F

(43)

in the cost function. Using control theory, the governing equations of the ﬂowﬁeld are introduced as a constraint in such a way that the ﬁnal expression for the gradient does not require reevaluation of the ﬂowﬁeld. In order to achieve this δw must be eliminated from (43). Suppose that the governing equation R which expresses the dependence of w and F within the ﬂowﬁeld domain D can be written as

In reference [72] the author derived the adjoint equations for transonic ﬂows modelled by both the potential ﬂow equation and the Euler equations. The theory was developed in terms of partial diﬀerential equations, leading to an adjoint partial diﬀerential equation. In order to obtain numerical solutions both the ﬂow and the adjoint equations must be discretized. The control theory might be applied directly to the discrete ﬂow equations which result R (w, F ) = 0. (44) from the numerical approximation of the ﬂow equations by ﬁnite element, ﬁnite volume or ﬁnite difThen δw is determined from the equation ference procedures. This leads directly to a set of discrete adjoint equations with a matrix which is ∂R ∂R δR = δw + δF = 0. (45) the transpose of the Jacobian matrix of the full set ∂w ∂F of discrete nonlinear ﬂow equations. On a threeNext, introducing a Lagrange Multiplier ψ, we have dimensional mesh with indices i, j, k the individual adjoint equations may be derived by collecting

∂I T ∂I T ∂R ∂R together all the terms multiplied by the variation δw + δF − ψ T δI = δw + δF ∂w ∂F ∂w ∂F T δwi,j,k of the discrete ﬂow variable wi,j,k . The re ∂I ∂R ∂I T ∂R = δw + δF. sulting discrete adjoint equations represent a possi− ψT − ψT ∂w ∂w ∂F ∂F ble discretization of the adjoint partial diﬀerential equation. If these equations are solved exactly they Choosing ψ to satisfy the adjoint equation can provide an exact gradient of the inexact cost T function which results from the discretization of the ∂R ∂I (46) ψ= ﬂow equations. On the other hand any consistent ∂w ∂w discretization of the adjoint partial diﬀerential equathe ﬁrst term is eliminated, and we ﬁnd that tion will yield the exact gradient in the limit as the mesh is reﬁned. The trade-oﬀ between the complexδI = GδF , (47) ity of the adjoint discretization, the accuracy of the resulting estimate of the gradient, and its impact where T on the computational cost to approach an optimum ∂I ∂R G= − ψT . solution is a subject of ongoing research, as indeed ∂F ∂F also remains true for the discretization of the ﬂow The advantage is that (47) is independent of δw, equations. with the result that the gradient of I with respect The true optimum shape belongs to an inﬁnitely dito an arbitrary number of design variables can be demensional space of design parameters. One motitermined without the need for additional ﬂow-ﬁeld vation for developing the theory for the partial difevaluations. In the case that (44) is a partial diﬀerferential equations of the ﬂow is to provide an inential equation, the adjoint equation (46) is also a dication in principle of how such a solution could partial diﬀerential equation and appropriate boundbe approached if suﬃcient computational resources ary conditions must be determined. were available. Another motivation is that it highAfter making a step in the negative gradient direclights the possibility of generating ill posed formulation, the gradient can be recalculated and the protions of the problem. For example, if one attempts cess repeated to follow a path of steepest descent to calculate the sensitivity of the pressure at a paruntil a minimum is reached. In order to avoid viticular location to changes in the boundary shape, olating constraints, such as a minimum acceptable there is the possibility that a shape modiﬁcation wing thickness, the gradient may be projected into could cause a shock wave to pass over that location. the allowable subspace within which the constraints Then the sensitivity could become unbounded. The are satisﬁed. In this way one can devise procedures movement of the shock, however, is continuous as which must necessarily converge at least to a local the shape changes. Therefore a quantity such as the minimum, and which can be accelerated by the use drag coeﬃcient, which is determined by integrating of more sophisticated descent methods such as conthe pressure over the surface, also depends continujugate gradient or quasi-Newton algorithms. There ously on the shape. The adjoint equation allows the is the possibility of more than one local minimum, sensitivity of the drag coeﬃcient to be determined but in any case the method will lead to an improvewithout the explicit evaluation of pressure sensitiviment over the original design. Furthermore, unlike ties which would be ill posed. the traditional inverse algorithms, any measure of The discrete adjoint equations, whether they are deperformance can be used as the cost function. 24

rived directly or by discretization of the adjoint partial diﬀerential equation, are linear. Therefore they could be solved by direct numerical inversion. In three-dimensional problems on a mesh with, say, n intervals in each coordinate direction, the number of unknowns is proportional to n3 and the bandwidth to n2 . The complexity of direct inversion is proportional to the number of unknowns multiplied by the square of the bandwidth, resulting in a complexity proportional to n7 . The cost of direct inversion can thus become prohibitive as the mesh is reﬁned, and it becomes more eﬃcient to use iterative solution methods. Moreover, because of the similarity of the adjoint equations to the ﬂow equations, the same iterative methods which have been proved to be eﬃcient for the solution of the ﬂow equations are eﬃcient for the solution of the adjoint equations.

the x1 , x2 , and x3 directions. Also introduce scaled contravariant velocity components

The control theory formulation for optimal aerodynamic design has proved eﬀective in a variety of applications [73, 77, 142]. The adjoint equations have also been used by Ta’asan, Kuruvila and Salas [166], who have implemented a one shot approach in which the constraint represented by the ﬂow equations is only required to be satisﬁed by the ﬁnal converged solution, and computational costs are also reduced by applying multigrid techniques to the geometry modiﬁcations as well as the solution of the ﬂow and adjoint equations. Pironneau has studied the use of control theory for optimal shape design of systems governed by elliptic equations [134], and more recently the Navier-Stokes equations, and also wave reﬂection problems. Adjoint methods have also been used by Baysal and Eleshaky [16].

Assume now that the new computational coordinate system conforms to the wing in such a way that the wing surface BW is represented by ξ2 = 0. Then the ﬂow is determined as the steady state solution of equation (48) subject to the ﬂow tangency condition

6.3

Ui = Qij uj . The transformed equations can now be written as ∂Fi ∂W + =0 ∂t ∂ξi

(48)

where W = Jw and

Fi = Qij fj =

ρUi ρUi u1 + Qi1 p ρUi u2 + Qi2 p ρUi u3 + Qi3 p ρUi H

U2 = 0 on BW .

(49)

At the far ﬁeld boundary BF , conditions are speciﬁed for incoming waves, as in the two-dimensional case, while outgoing waves are determined by the solution. The weak form of the Euler equations for steady ﬂow can be written as ∂φT Fi dD = ni φT Fi dB, (50) D ∂ξi B where the test vector φ is an arbitrary diﬀerentiable function and ni is the outward normal at the boundary. If a diﬀerentiable solution w is obtained to this equation, it can be integrated by parts to give ∂Fi φT dD = 0 ∂ξi D

Three-Dimensional Design using the Euler Equations

In order to illustrate the application of control theory to aerodynamic design problems, this section treats the case of three-dimensional wing design using the inviscid Euler equations as the mathematical model for compressible ﬂow. A transformation to a body-ﬁtted coordinate system will be introduced, so that variations in the wing shape induce corresponding variations in the computational mesh. Thus the ﬂow is determined by the solution of the transformed equation (5). Let ∂xi ∂ξi −1 = , J = det(K), Kij , Kij = ∂ξj ∂xj

and since this is true for any φ, the diﬀerential form can be recovered. If the solution is discontinuous, equation (50) may be integrated by parts separately on either side of the discontinuity to recover the shock jump conditions. Suppose now that it is desired to control the surface pressure by varying the wing shape. It is convenient to retain a ﬁxed computational domain. Variations in the shape then result in corresponding variations in the mapping derivatives deﬁned by K. Introduce the cost function 1 2 (p − pd ) dξ1 dξ3 , I= 2 BW

and Q = JK −1 . The elements of Q are the coeﬃcients of K, and in a ﬁnite volume discretization they are just the face areas of the computational cells projected in

where pd is the desired pressure. The design problem is now treated as a control problem where the control 25

function is the wing shape, which is to be chosen to minimize I subject to the constraints deﬁned by the ﬂow equations (48–50). A variation in the shape will cause a variation δp in the pressure and consequently a variation in the cost function (p − pd ) δp dξ1 dξ3 . (51) δI =

Then if the coordinate transformation is such that δQ is negligible in the far ﬁeld, the only remaining boundary term is ψ T δF2 dξ1 dξ3 . − BW

Thus by letting ψ satisfy the boundary condition,

BW

Q21 ψ2 + Q22 ψ3 + Q23 ψ4 = (p − pd ) on BW ,

Since p depends on w through the equation of state (2), the variation δp can be determined from the variation δw. Deﬁne the Jacobian matrices Ai =

∂fi , ∂w

we ﬁnd ﬁnally that δI

Ci = Qij Aj .

(56)

(52)

=− D

−

The weak form of the equation for δw in the steady state becomes ∂φT δFi dD = (ni φT δFi )dB, D ∂ξi B

∂ψ T δQij fj dD ∂ξi

(δQ21 ψ2 + δQ22 ψ3 + Q23 ψ4 ) p dξ1 dξ3 . (57) BW

A convenient way to treat a wing is to introduce η

where δFi = Ci δw + δQij fj ,

s ξ

which should hold for any diﬀerential test function φ. This equation may be added to the variation in the cost function, which may now be written as δI = (p − pd ) δp dξ1 dξ3

15a: x, y-Plane.

Figure 15: Sheared Parabolic Mapping. sheared parabolic coordinates as shown in ﬁgure 15 through the transformation 1 x = x0 (ζ) + a (ζ) ξ 2 − (η + S (ξ, ζ))2 2 y = y0 (ζ) + a (ζ) ξ (η + S (ξ, ζ)) z = ζ.

BW

∂ψ T δFi dD ∂ξi D ni ψ T δFi dB. +

−

(53)

B

On the wing surface BW , from equation (49) that 0 Q21 δp δF2 = Q22 δp Q23 δp 0

n1 = n3 = 0 and it follows

Here x = x1 , y = x2 , z = x3 are the Cartesian coordinates, and ξ and η +S correspond to parabolic coordinates generated by the mapping

0

δQ21 p + δQ22 p . δQ23 p 0

1 x + iy = x0 + iy0 + a (ζ) {ξ + i (η + S)}2 2

(54)

at a ﬁxed span station ζ. x0 (ζ) and y0 (ζ) are the coordinates of a singular line which is swept to lie just inside the leading edge of a swept wing, while a (ζ) is a scale factor to allow for spanwise chord variations.

Since the weak equation for δw should hold for an arbitrary choice of the test vector φ, we are free to choose φ to simplify the resulting expressions. Therefore we set φ = ψ, where the costate vector ψ is the solution of the adjoint equation ∂ψ ∂ψ − CiT = 0 in D. ∂t ∂ξi

15b: ξ, η-Plane.

We now treat S (ξ, ζ) as the control. Substitution of these formulas yields the variation in the form δI = G(ξ, η)δS(ξ, η)dξdη

(55)

where the gradient G(ξ, η) is obtained by evaluating the integrals in equation (57). Thus to reduce I we can choose δS = −λG

At the outer boundary incoming characteristics for ψ correspond to outgoing characteristics for δw. Consequently one can choose boundary conditions for ψ such that ni ψ T Ci δw = 0.

where λ is suﬃciently small and non-negative. In order to impose a thickness constraint we can deﬁne 26

a baseline surface S0 (ξ, ζ) below which S (ξ, ζ) is not allowed to fall. Now we take λ = λ (ξ, ζ) as a non-negative function such that

Thus the gradient with respect to the B-spline coefﬁcients is obtained by multiplying G by B T , and a descent step is deﬁned by setting

S (ξ, ζ) + δS (ξ, ζ) ≥ S0 (ξ, ζ) .

d = −λB T G, δS = Bd = −λBB T G

(58)

Then the constraint is satisﬁed, while δI = − λG 2 dξ dζ ≤ 0.

where λ is suﬃciently small and positive. The coeﬃcients of B can be renormalized to produce unit row sums. With a uniform mesh spacing in the computational domain this formula is equivalent to the use of a gradient modiﬁed by two passes of the explicit smoothing procedure

BW

The costate solution ψ is a legitimate test function for the weak form of the ﬂow equations only if it is diﬀerentiable. Smoothness should also be preserved in the redesigned shape. It is therefore crucially important to introduce appropriate smoothing procedures. In order to avoid discontinuities in the adjoint boundary condition which would be caused by the appearance of shock waves, the cost function for the target pressure may be modiﬁed to the form 2 ' & ∂Z 1 dξdη λ1 Z + λ2 I = 2 ∂ξ λ1 Z Then δI

−

1 2 1 G¯i,k = Gi−1,k + Gi,k + Gi+1,k 6 3 6 with a similar smoothing procedure in the k discretization. Implicit smoothing may also be used. The smoothing equation −i+ 12 ,k (G¯i+1,k − G¯i,k ) + i− 12 ,k (G¯i,k − G¯i−1,k ) = Gi,k approximates the diﬀerential equation

∂ ∂Z λ2 = p − pd . ∂ξ ∂ξ

∂ ∂ G¯ =G G¯ − ∂ξ ∂ξ

∂Z ∂ δZ dξdη = λ1 ZδZ + λ2 ∂ξ ∂ξ ∂ ∂ λ2 = Z λ1 − δZ dξdη ∂ξ ∂ξ = Zδp dξdη

¯ then to ﬁrst order the change If one sets δS = −λG, in the cost is δI = − GδS dξdη 1 ∂ ∂ G¯ ¯ = − G¯ − G dξdη λ ∂ξ ∂ξ ¯ 2 ' & ∂G 1 2 G¯ + = − dξdη λ ∂ξ

and the smooth quantity Z replaces p − pd in the adjoint boundary condition.

< 0,

Independent movement of the boundary mesh points could produce discontinuities in the designed shape. In order to prevent this the gradient may be also smoothed. Both explicit and implicit smoothing procedures are useful. Suppose that the movement of the surface mesh points were deﬁned by local Bsplines. In the case of a uniform one-dimensional mesh, a B-spline with a displacement d centered at the mesh point i would produce displacements d/4 at i+1 and i−1 and zero elsewhere, while preserving continuity of the ﬁrst and second derivatives. Thus we can suppose that the discrete surface displacement has the form

assuring an improvement if λ is suﬃciently small and positive, unless the process has already reached a stationary point at which G = 0.

6.4

Design of Swept Wings for Very Low Shock Drag

The method has been used to carry out a study of swept wing designs which might be appropriate for long range transport aircraft. Since three dimensional calculations require substantial computational resources, it is extremely important for the practical implementation of the method to use fast solution algorithms for the ﬂow and the adjoint equations. In this case the author’s FLO87 computer program has been used as the basis of the design method. FLO87 solves the three dimensional Euler equations with a cell-centered ﬁnite volume scheme, and uses residual averaging and multigrid acceleration to obtain

δS = Bd, where B is a matrix with coeﬃcients deﬁned by the B-splines, and di is the displacement associated with the B-spline centered at i. Then, using the discrete formulas, to ﬁrst order the change in the cost is δI = G T δS = G T Bd. 27

very rapid steady state solutions, usually in 25 to 50 multigrid cycles [66, 70]. Upwind biasing is used to produce non-oscillatory solutions, and assure the clean capture of shock waves. This is introduced through the addition of carefully controlled numerical diﬀusion terms, with a magnitude of order ∆x3 in smooth parts of the ﬂow. The adjoint equations are treated in the same way as the ﬂow equations. The ﬂuxes are ﬁrst estimated by central diﬀerences, and then modiﬁed by downwind biasing through numerical diﬀusive terms which are supplied by the same subroutines that were used for the ﬂow equations.

ﬁcient to provide a shape modiﬁcation which leads to an improvement. Figures 23 and 24 show a wing which was designed for a lift coeﬃcient of .50 at Mach .85. In order to prevent the ﬁnal wing from becoming too thin the threshold S0 (ξ, η) was set at three quarters of the height of the bump S (ξ, η) deﬁning the initial wing. This calculation was performed on a mesh with 192 intervals in the ξ direction wrapping around the wing, 32 intervals in the normal η direction and 48 intervals in the spanwise ζ direction, giving a total of 294912 cells. The wing was speciﬁed by 33 sections, each with 128 points, giving a total of 4224 design variables. The plots show the initial wing geometry and pressure distribution, and the modiﬁed geometry and pressure distribution after 40 design cycles. The total inviscid drag coeﬃcient was reduced from 0.0210 to 0.0112. The initial design exhibits a very strong shock wave in the inboard region. It can be seen that this is completely eliminated, leaving a very weak shock wave in the outboard region. To verify the solution, the ﬁnal geometry was analyzed with another method, using the computer program FLO67. This program uses a cell-vertex formulation, and has recently been modiﬁed to incorporate a local extremum diminishing algorithm with a very low level of numerical diﬀusion [76]. When run to full convergence it was found that a better estimate of the drag coeﬃcient of the redesigned wing is 0.0094 at Mach 0.85 with a lift coeﬃcient of 0.5, giving a lift to drag ratio of 53. The results from FLO67 for the initial and ﬁnal wings are illustrated in Figures 25 and 26. A calculation at Mach 0.500 shows a drag coeﬃcient of 0.0087 for a lift coeﬃcient of 0.5. Since in this case the ﬂow 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 drag, CD = CL 2 /πAR, with an aspect ratio AR = 9, and an eﬃciency factor = 0.97. Thus the design method has reduced the shock wave drag coeﬃcient to about 0.0007 at a lift coeﬃcient of 0.5. Figure 27 shows the result of an analysis for an oﬀ design point with the Mach number increased to .86 with the same lift coeﬃcient of .5. This results in a ﬂat-topped pressure distribution terminating with a weak shock of nearly uniform strength across the whole span. The drag coeﬃcient is .0097. The penalty of .0003 is so small that this might be a preferred cruising condition.

The study has been focussed on wings designed for cruising at Mach .85, with lift coeﬃcients in the range of .5 to .55. In every case, the wing planform was ﬁxed while the sections were free to be changed arbitrarily by the design method, with a restriction on the minimum thickness. The initial wing has a unit-semi-span, with 38 degrees leading edge sweep. It has a modiﬁed trapezoidal planform, with straight taper from a root chord of 0.38, and a curved trailing edge in the inboard region blending into straight taper outboard of the 30 percent span station to a tip chord of 0.10, with an aspect ratio of 9.0. The initial wing sections were based on a section speciﬁcally designed by the author’s two dimensional design method [73] to give shock free ﬂow at Mach 0.78 with a lift coeﬃcient of 0.6. This section, which has a thickness to chord ratio of 9.5 percent, was used at the tip. Similar sections with an increased thickness were used inboard. The variation of thickness was non-linear with a more rapid increase near the root, where the thickness to chord ratio 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 from root to tip. The two-dimensional pressure distribution of the basic wing section at its design point was introduced as a target pressure distribution uniformly across the span. This target is presumably not realizable, but serves to favor the establishment of a relatively benign pressure distribution. The total inviscid drag coeﬃcient, due to the combination of vortex and shock wave drag, was also included in the cost function. Since the main objective of the study was to minimize the drag, the target pressure distribution was reset after every fourth design cycle to a distribution derived by smoothing the existing pressure distribution. This allows the scheme more freedom to make changes which reduce drag. The calculations were performed with the lift coeﬃcient forced to approach a ﬁxed value by adjusting the angle of attack every ﬁfth iteration of the ﬂow solution. It was found that the computational costs can be reduced by using only 15 multigrid cycles in each ﬂow solution, and in each adjoint solution. Although this is not enough for full convergence, it proves suf-

A second wing was designed in exactly the same manner as the ﬁrst, starting from the same initial geometry and with the same constraints, to give a lift coeﬃcient of .55 at Mach .85. This produces stronger shock waves and is therefore a more severe test of the method. In this case the total inviscid drag coeﬃcient was reduced from 0.0243 to 0.0134 28

in 40 design cycles. Again the performance of the ﬁnal design was veriﬁed by a calculation with FLO67, and when the result was fully converged the drag coeﬃcient was found to be 0.0115. A subsonic calculation at Mach .500 shows a drag coeﬃcient of 0.0107 for a lift coeﬃcient of 0.55. Thus in this case the shock wave drag coeﬃcient is about 0.0008. For a representative transport aircraft the parasite drag coeﬃcient of the wing due to skin friction is about 0.0045. Also the fuselage drag coeﬃcient is about 0.0050, the nacelle drag coeﬃcient is about 0.0015, the empennage drag coeﬃcient is about 0.0020, and excrescence drag coeﬃcient is about 0.0010. This would give a total drag coeﬃcient CD = 0.0255 for a lift coeﬃcient of 0.55, corresponding to a lift to drag ratio L/D = 21.6. This would be a substantial improvement over the values obtained by currently ﬂying transport aircraft.

A detailed derivation is given in reference [78]. Thus the perturbation equation can be written as ∂fj ∂ (δw + δw∗ ) = 0 ∂ξi ∂w where δw is the variation in the solution at ﬁxed ξ caused by the change in the boundary, while δw∗ is the change in the original solution w(ξ) corresponding to the mesh movement δx(ξ)

Now ψT

∂ δFi dξ ∂ξi

= B

−

6.5

Optimization of Complex Conﬁgurations

= R(η)δxb (ξ)

δy(ξ, η)

= R(η)δyb (ξ)

D

ni ψ T δFi dξB ∂ψ T Ci (δw + δw∗ )dξ ∂ξi

and if ψ satisﬁes the adjoint equation the entire ﬁeld integral is eliminated, leaving only the boundary integral in equation (57).

In order to treat more complex conﬁgurations one can use a numerical grid generation procedure to produce a body-ﬁtted mesh for the initial geometry, and then modify the mesh in subsequent design cycles by an analytic perturbation formula. In the twodimensional case, for example, with computational coordinates ξ, η, let the boundary displacement at η = 0 be δxb (ξ), δyb (ξ). Then the mesh points along the radial coordinate lines ξ = constant can be replaced by δx(ξ, η)

∂wi δξk ∂ξk

δwi∗ =

In an actual discretization the ﬁeld terms are not zero, but this result suggests that they should be small is a ﬁne enough mesh is used, and might be dropped. This allows a drastic simpliﬁcation of the treatment of complex conﬁgurations. Preliminary numerical experiments with airfoil and wing calculations indicate roughly the same convergence with and without the ﬁeld terms in the gradient.

7

OUTLOOK CONCLUSIONS

AND

yielding ( δK =

∂ R(η) ∂ξ δxb ∂ R(η) ∂ξ δyb

dR dη δxb dR dη δyb

)

Better algorithms and better computer hardware have contributed about equally to the progress of computational science in the last two decades. In 1970 the Control Data 6600 represented the state of the art in computer hardware with a speed of about 106 operations per second (one megaﬂop), while in 1990 the 8 processor Cray YMP oﬀered a performance of about 109 operations per second (one gigaﬂop). Correspondingly, steady-state Euler calculations which required 5,000–10,000 steps prior to 1980 could be performed in 10–50 steps in 1990 using multigrid acceleration. With the advent of massively parallel computers it appears that the progress of computer hardware may even accelerate. Teraﬂop machines oﬀering further improvement by a factor of 1,000 are likely to be available within a few years. Parallel architectures will force a reappraisal of existing algorithms, and their eﬀective utilization will require the extensive development of new parallel software.

Such a procedure has been implemented by J. Reuther for the three-dimensional Euler equations, and applied to the optimization of wing-body conﬁgurations [143]. It is also possible to show that in the continuous limit the ﬁeld integral in equation (57) can be eliminated. Let the change in the coordinates xi at ﬁxed ξ be δxi (ξ). Then, using the fact that the ﬂuxes fj (w) satisfy the ﬂow equation (48), it is possible to show by a direct calculation that ∂ ∂fj ∂w ∂ δQij fj = Qij δξk ∂ξi ∂ξi ∂w ∂ξk where

δξ = K −1 δx. 29

In parallel with the transition to more sophisticated algorithms, the present challenge is to extend the effective use of CFD to more complex applications. A key problem is the treatment of multiple space and time scales. These arise not only in turbulent ﬂows, but also in many other situations such as chemically reacting ﬂows, combustion, ﬂame fronts and plasma dynamics. Another challenge, is presented by problems with moving boundaries. Examples include helicopter rotors, and rotor-stator interaction in turbomachinery. Algorithms for these problems can be signiﬁcantly improved by innovative concepts, such as the idea of time inclining. It can be anticipated that interdisciplinary applications in which CFD is coupled with the computational analysis of other properties of the design will play an increasingly important role. These applications may include structural, thermal and electromagnetic analysis. Aeroelastic problems and integrated control system and aerodynamic design are likely target areas. The

Initial Design

CAD Definition

CFD Geometry Modeling

the way in which a numerical wind tunnel might evolve from current techniques, which involve massive data handling tasks, to a fully integrated design environment. In the long run, computational simulation should become the principal tool of the aerodynamic design process because of the ﬂexibility it provides for the rapid and comparatively inexpensive evaluation of alternative designs, and because it can be integrated in a numerical design environment providing for both multi-disciplinary analysis and multidisciplinary optimization.

REFERENCES [1] R. Abid, C.G. Speziale, and S. Thangam. Application of a new k-τ model to near wall turbulent ﬂows. AIAA Paper 91-0614, AIAA 29th Aerospace Sciences Meeting, Reno, NV, January 1991.

Flow Solution

Mesh Generation

[2] H. Aiso. Admissibility of diﬀerence approximations for scalar conservation laws. Hiroshima Math. Journal, 23:15–61, 1993.

Aeroelastic Solution

Requirements

Redesign

Quantitative Assessment L, D, W MDO*

Loads

Visualization

[3] J. J. Alonso and A. Jameson. Fully-implicit time-marching aeroelastic solutions. AIAA paper 94-0056, AIAA 32nd Aerospace Sciences Meeting, Reno, Nevada, January 1994.

Numerically Intensive

Human Intensive

*MDO: Multi-Disciplinary Optimization

[4] J. J. Alonso, L. Martinelli, and A. Jameson. Multigrid unsteady Navier-Stokes calculations with aeroelastic applications. AIAA paper 950048, AIAA 33rd Aerospace Sciences Meeting, Reno, Nevada, January 1995.

Figure 16: Concept for a numerical wind tunnel.

Initial Design

Master Definition Central Database

Automatic Mesh Generation

Geometry Modification

[5] B.K. Anderson, J.L. Thomas, and B. Van Leer. A comparison of ﬂux vector splittings for the Euler equations. AIAA Paper 85-0122, Reno, NV, January 1985.

Flow Solution Optimization Aeroelastic Solution Requirements

High-Level Redesign

Human Intensive

Monitor Results

Loads

Quantitative Assessment

Numerically Intensive

[6] W.K. Anderson, J.L. Thomas, and D.L. Whitﬁeld. Multigrid acceleration of the ﬂux split Euler equations. AIAA Paper 86-0274, AIAA 24th Aerospace Sciences Meeting, Reno, January 1986.

Figure 17: Advanced numerical wind tunnel. development of improved algorithms continues to be important in providing the basic building blocks for numerical simulation. In particular, better error estimation procedures must be developed and incorporated in the simulation software to provide error control. The basic simulation software is only one of the needed ingredients, however. The ﬂow solver must be embedded in a user-friendly system for geometry modeling, output analysis, and data management that will provide a complete numerical design environment. These are the ingredients which are needed for the full realization of the concept of a numerical wind tunnel. Figures 16 and 17 illustrate

[7] T.J. Baker. Mesh generation by a sequence of transformations. Appl. Num. Math., 2:515– 528, 1986. [8] N. Balakrishnan and S. M. Deshpande. New upwind schemes with wave-particle splitting for inviscid compressible ﬂows. Report 91 FM 12, Indian Institute of Science, 1991. [9] B. Baldwin and H. Lomax. Thin layer approximation and algebraic model for separated turbulent ﬂow. AIAA Paper 78-257, 1978. 30

[10] B.S. Baldwin and T.J. Barth. A one-equation turbulence transport model for high Reynolds number wall-bounded ﬂows. AIAA Paper 910610, AIAA 29th Aerospace Sciences Meeting, Reno, NV, January 1991.

[21] M. Berger and A. Jameson. Automatic adaptive grid reﬁnement for the Euler equations. AIAA Journal, 23:561–568, 1985. [22] M. Berger and R.J. LeVeque. An adaptive Cartesian mesh algorithm for the Euler equations in arbitrary geometries. AIAA Paper 891930, 1989.

[11] T. J. Barth. Aspects of unstructured grids and ﬁnite volume solvers for the Euler and Navier Stokes equations. In von Karman Institute for Fluid Dynamics Lecture Series Notes 1994-05, Brussels, 1994.

[23] J.P. Boris and D.L. Book. Flux corrected transport, 1 SHASTA, a ﬂuid transport algorithm that works. J. Comp. Phys., 11:38–69, 1973.

[12] T.J. Barth and P.O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, January 1990.

[24] A. Brandt. Multi-level adaptive solutions to boundary value problems. Math. Comp., 31:333–390, 1977. [25] M.O. Bristeau, R. Glowinski, J. Periaux, P. Perrier, O. Pironneau, and C. Poirier. On the numerical solution of nonlinear problems in ﬂuid dynamics by least squares and ﬁnite element methods (II), application to transonic ﬂow simulations. Comp. Meth. Appl. Mech. and Eng., 51:363–394, 1985.

[13] T.J. Barth and D.C. Jespersen. The design and application of upwind schemes on unstructured meshes. AIAA paper 89-0366, AIAA 27th Aerospace Sciences Meeting, Reno, Nevada, January 1989. [14] J.T. Batina. Implicit ﬂux-split Euler schemes for unsteady aerodynamic analysis involving unstructured dynamic meshes. AIAA paper 90-0936, April 1990.

[26] R.J. Busch, Jr. Computational ﬂuid dynamics in the design of the Northrop/McDonnell Douglas YF-23 ATF prototype. AIAA paper 91-1627, AIAA 21st Fluid Dynamics, Plasmadynamics & Lasers Conference, Honolulu, Hawaii, 1991.

[15] F. Bauer, P. Garabedian, D. Korn, and A. Jameson. Supercritical Wing Sections II. Springer Verlag, New York, 1975.

[27] C. Canuto, M.Y. Hussaini, A. Quarteroni, and D.A. Zang. Spectral Methods in Fluid Dynamics. Springer Verlag, 1987.

[16] O. Baysal and M. E. Eleshaky. Aerodynamic design optimization using sensitivity anaysis and computational ﬂuid dynamics. AIAA paper 91-0471, 29th Aerospace Sciences Meeting, Reno, Nevada, January 1991.

[28] M.H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable boundary conditions for ﬁnite-diﬀerence schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. ICASE Report 93-9, Hampton, VA, March 1993.

[17] R.W. Beam and R.F. Warming. An implicit ﬁnite diﬀerence algorithm for hyperbolic systems in conservation form. J. Comp. Phys., 23:87–110, 1976.

[29] D.A. Caughey. A diagonal implicit multigrid algorithm for the Euler equations. AIAA Paper 87-453, 25th Aerospace Sciences Meeting, Reno, January 1987.

[18] A. Belov, L. Martinelli, and A. Jameson. A new implicit algorithm with multigrid for unsteady incompressible ﬂow calculations. AIAA paper 95-0049, AIAA 33rd Aerospace Sciences Meeting, Reno, Nevada, January 1995.

[30] T. Cebeci and A.M.O. Smith. Analysis of Turbulent Boundary Layers. Academic Press, 1974.

[19] J.A. Benek, P.G. Buning, and J.L. Steger. A 3-D Chimera grid embedding technique. In Proceedings AIAA 7th Computational Fluid Dynamics Conference, pages 507–512, Cincinnati, OH, 1985. AIAA Paper 85-1523.

[31] S.R. Chakravarthy. Relaxation methods for unfactored implicit upwind schemes. AIAA Paper 84-0165, AIAA 22nd Aerospace Sciences Meeting, Reno, January 1984.

[20] J.A. Benek, T.L. Donegan, and N.E. Suhs. Extended Chimera grid embedding scheme with applications to viscous ﬂows. AIAA Paper 871126, AIAA 8th Computational Fluid Dynamics Conference, Honolulu, HI, 1987.

[32] S.R. Chakravarthy, A. Harten, and S. Osher. Essentially non-oscillatory shock capturing schemes of uniformly very high accuracy. AIAA Paper 86-0339, AIAA 24th Aerospace Sciences Meeting, Reno, January 1986. 31

[33] R. Chipman and A. Jameson. Fully conservative numerical solutions for unsteady irrotational transonic ﬂow about airfoils. AIAA Paper 79-1555, AIAA 12th Fluid and Plasma Dynamics Conference, Williamsburg, VA, July 1979.

[45] W. Eppard and B. Grossman. A multidimensional kinetic-based upwind solver for the Euler equations. AIAA Paper 93-3303, AIAA 11th Computational Fluid Dynamics Conference, Orlando, FL, July 1993. [46] L.E. Eriksson. Generation of boundaryconforming grids around wing-body conﬁgurations using transﬁnite interpolation. AIAA Journal, 20:1313–1320, 1982.

[34] S.E. Cliﬀ and S.D. Thomas. Euler/experimental correlations of sonic boom pressure signatures. AIAA Paper 91-3276, AIAA 9th Applied Aerodynamics Conference, Baltimore, September 1991.

[47] J. Farmer, L. Martinelli, A. Jameson, and G. Cowles. Fully-nonlinear CFD techniques for ship performance analysis and design. AIAA paper 95-1690, AIAA 12th Computational Fluid Dynamics Conference, San Diego, CA, June 1995.

[35] T.J. Coakley. Numerical simulation of viscous transonic airfoil ﬂows. AIAA Paper 87-0416, AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.

[48] R.P. Fedorenko. The speed of convergence of one iterative process. USSR Comp. Math. and Math. Phys., 4:227–235, 1964.

[36] J.P. Croisille and P. Villedieu. Kinetic ﬂux splitting schemes for hypersonic ﬂows. In M. Napolitano and F. Sobetta, editors, Proc 13th International Congress on Numerical Methods in Fluid Dynamics, pages 310–313, Rome, July 1992. Springer Verlag.

[49] C.W. Gear. The numerical integration of stiﬀ ordinary diﬀerential equations. Report 221, University of Illinois Department of Computer Science, 1967.

[37] R.M. Cummings, Y.M. Rizk, L.B. Schiﬀ, and N.M. Chaderjian. Navier-Stokes predictions for the F-18 wing and fuselage at large incidence. J. of Aircraft, 29:565–574, 1992.

[50] M. Giles, M. Drela, and W.T. Thompkins. Newton solution of direct and inverse transonic Euler equations. AIAA Paper 85-1530, Cincinnati, 1985.

[38] G. Dahlquist. A special stability problem for linear multistep methods. BIT, 3:27–43, 1963.

[51] S.K. Godunov. A diﬀerence method for the numerical calculation of discontinuous solutions of hydrodynamic equations. Mat. Sbornik, 47:271–306, 1959. Translated as JPRS 7225 by U.S. Dept. of Commerce, 1960.

[39] J.F. Dannenhoﬀer and J.R. Baron. Robust grid adaptation for complex transonic ﬂows. AIAA Paper 86-0495, AIAA 24th Aerospace Sciences Meeting, Reno, January 1986.

[52] M.H. Ha. The impact of turbulence modelling on the numerical prediction of ﬂows. In M. Napolitano and F. Solbetta, editors, Proc. of the 13th International Conference on Numerical Methods in Fluid Dynamics, pages 27– 46, Rome, Italy, July 1992. Springer Verlag, 1993.

[40] D. Deganiand and L. Schiﬀ. Computation of turbulent supersonic ﬂows around pointed bodies having crossﬂow separation. J. Comp. Phys., 66:173–196, 1986. [41] B. Delaunay. Sur la sph`ere vide. Bull. Acad. Science USSR VII: Class Scil, Mat. Nat., pages 793–800, 1934.

[53] W. Hackbusch. On the multi-grid method applied to diﬀerence equations. Computing, 20:291–306, 1978.

[42] S.M. Deshpande. On the Maxwellian distribution, symmetric form and entropy conservation for the Euler equations. NASA TP 2583, 1986.

[54] M. Hafez, J.C. South, and E.M. Murman. Artiﬁcial compressibility method for numerical solutions of the transonic full potential equation. AIAA Journal, 17:838–844, 1979.

[43] A. Eberle. A ﬁnite volume method for calculating transonic potential ﬂow around wings from the minimum pressure integral. NASA TM 75324, 1978. Translated from MBB UFE 1407(0).

[55] M.G. Hall. Cell vertex multigrid schemes for solution of the Euler equations. In Proc. IMA Conference on Numerical Methods for Fluid Dynamics, Reading, April 1985.

[44] P.R. Eiseman. A multi-surface method of coordinate generation. J. Comp. Phys., 33:118– 150, 1979.

[56] A. Harten. High resolution schemes for hyperbolic conservation laws. J. Comp. Phys., 49:357–393, 1983. 32

[57] P.W. Hemker and S.P. Spekreijse. Multigrid solution of the steady Euler equations. In Proc. Oberwolfach Meeting on Multigrid Methods, December 1984.

October 1985. Princeton University Report MAE 1743. [68] A. Jameson. Non-oscillatory shock capturing scheme using ﬂux limited dissipation. In B.E. Engquist, S. Osher, and R.C.J. Sommerville, editors, Lectures in Applied Mathematics, Vol. 22, Part 1, Large Scale Computations in Fluid Mechanics, pages 345–370. AMS, 1985.

[58] J.L. Hess and A.M.O. Smith. Calculation of non-lifting potential ﬂow about arbitrary three-dimensional bodies. Douglas Aircraft Report ES 40622, 1962. [59] C. Hirsch, C. Lacol, and H. Deconinck. Convection algorithms based on a diagonalization procedure for the multi-dimensional Euler equations. In Proc. AIAA 8th Computational Fluid Dynamics Conference, pages 667– 676, Hawaii, June 1987. AIAA Paper 87-1163.

[69] A. Jameson. Transonic ﬂow calculations for aircraft. In F. Brezzi, editor, Lecture Notes in Mathematics, Numerical Methods in Fluid Dynamics, pages 156–242. Springer Verlag, 1985. [70] A. Jameson. Multigrid algorithms for compressible ﬂow calculations. In W. Hackbusch and U. Trottenberg, editors, Lecture Notes in Mathematics, Vol. 1228, pages 166–201. Proceedings of the 2nd European Conference on Multigrid Methods, Cologne, 1985, SpringerVerlag, 1986.

[60] C. Hirsch and P. Van Ransbeeck. Multidimensional upwinding and artiﬁcial dissipation. Technical report. Published in Frontiers of Computational Fluid Dynamics 1994, D. A. Caughey and M. M. Hafez, editors, Wiley, pp. 597–626.

[71] A. Jameson. A vertex based multigrid algorithm for three-dimensional compressible ﬂow calculations. In T.E. Tezduar and T.J.R. Hughes, editors, Numerical Methods for Compressible Flow - Finite Diﬀerence, Element And Volume Techniques, 1986. ASME Publication AMD 78.

[61] D.G. Holmes and S.H. Lamson. Adaptive triangular meshes for compressible ﬂow solutions. In Proceedings First International Conference on Numerical Grid Generation in Computational Fluid Dynamics, pages 413– 424, Landshut, FRG, July 1986. [62] T.J.R. Hughes, L.P. Franca, and M. Mallet. A new ﬁnite element formulation for computational ﬂuid dynamics,I, Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics. Comp. Meth. Appl. Mech. and Eng., 59:223–231, 1986.

[72] A. Jameson. Aerodynamic design via control theory. J. Sci. Comp., 3:233–260, 1988. [73] A. Jameson. Automatic design of transonic airfoils to reduce the shock induced pressure drag. In Proceedings of the 31st Israel Annual Conference on Aviation and Aeronautics, Tel Aviv, pages 5–17, February 1990.

[63] A. Jameson. Iterative solution of transonic ﬂows over airfoils and wings, including ﬂows at Mach 1. Comm. on Pure and Appl. Math., 27:283–309, 1974.

[74] A. Jameson. Time dependent calculations using multigrid, with applications to unsteady ﬂows past airfoils and wings. AIAA paper 911596, AIAA 10th Computational Fluid Dynamics Conference, Honolulu, Hawaii, June 1991.

[64] A. Jameson. Transonic potential ﬂow calculations in conservation form. In Proc. AIAA 2nd Computational Fluid Dynamics Conference, pages 148–161, Hartford, 1975.

[75] A. Jameson. Artiﬁcial diﬀusion, upwind biasing, limiters and their eﬀect on accuracy and multigrid convergence in transonic and hypersonic ﬂows. AIAA Paper 93-3359, AIAA 11th Computational Fluid Dynamics Conference, Orlando, FL, July 1993.

[65] A. Jameson. Solution of the Euler equations by a multigrid method. Appl. Math. and Comp., 13:327–356, 1983. [66] A. Jameson. Solution of the Euler equations for two dimensional transonic ﬂow by a multigrid method. Appl. Math. Comp., 13:327–356, 1983.

[76] A. Jameson. Artiﬁcial diﬀusion, upwind biasing, limiters and their eﬀect on accuracy and multigrid convergence in transonic and hypersonic ﬂow. AIAA paper 93-3359, AIAA 11th Computational Fluid Dynamics Conference, Orlando, Florida, July 1993.

[67] A. Jameson. Multigrid algorithms for compressible ﬂow calculations. In Second European Conference on Multigrid Methods, Cologne, 33

[77] A. Jameson. Optimum aerodynamic design via boundary control. Technical report, AGARD FDP/Von Karman Institute Special Course on Optimum Design Methods in Aerodynamics, Brussels, April 1994.

[88] D. Johnson and L. King. A mathematically simple turbulence closure model for attached and separated turbulent boundary layers. AIAA Journal, 23:1684–1692, 1985. [89] W.P. Jones and B.E. Launder. The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence. Int. J. of Heat Tran., 16:1119–1130, 1973.

[78] A. Jameson. MAE Technical Report 2050, Princeton University, Princeton, New Jersey, October 1995.

[90] W.H. Jou. Boeing Memorandum AEROB113B-L92-018, September 1992. To Joseph Shang.

[79] A. Jameson. Analysis and design of numerical schemes for gas dynamics 1, artiﬁcial diﬀusion, upwind biasing, limiters and their eﬀect on multigrid convergence. Int. J. of Comp. Fluid Dyn., 4:171–218, 1995.

[91] T.J. Kao, T.Y. Su, and N.J. Yu. Navier-Stokes calculations for transport wing-body conﬁgurations with nacelles and struts. AIAA Paper 93-2945, AIAA 24th Fluid Dynamics Conference, Orlando, July 1993.

[80] A. Jameson. Analysis and design of numerical schemes for gas dynamics 2, artiﬁcial diﬀusion and discrete shock structure. Int. J. of Comp. Fluid Dyn., To Appear.

[92] G.E. Karniadakis and S.A. Orszag. Nodes, modes and ﬂow codes. Physics Today, pages 34–42, March 1993.

[81] A. Jameson and T.J. Baker. Improvements to the aircraft Euler method. AIAA Paper 870452, AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.

[93] M.H. Lallemand and A. Dervieux. A multigrid ﬁnite-element method for solving the twodimensional Euler equations. In S.F. McCormick, editor, Proceedings of the Third Copper Mountain Conference on Multigrid Methods, Lecture Notes in Pure and Applied Mathematics, pages 337–363, Copper Mountain, April 1987.

[82] A. Jameson, T.J. Baker, and N.P. Weatherill. Calculation of inviscid transonic ﬂow over a complete aircraft. AIAA Paper 86-0103, AIAA 24th Aerospace Sciences Meeting, Reno, January 1986. [83] A. Jameson and D.J. Mavriplis. Multigrid solution of the Euler equations on unstructured and adaptive grids. In S. McCormick, editor, Multigrid Methods, Theory, Applications and Supercomputing. Lecture Notes in Pure and Applied Mathematics, volume 110, pages 413– 430, April 1987.

[94] A.M. Landsberg, J.P. Boris, W. Sandberg, and T.R. Young. Naval ship superstructure design: Complex three-dimensional ﬂows using an efﬁcient, parallel method. High Performance Computing 1993: Grand Challenges in Computer Simulation, 1993. [95] P. D. Lax. Hyperbolic systems of conservation laws. SIAM Regional Series on Appl. Math., II, 1973.

[84] A. Jameson, W. Schmidt, and E. Turkel. Numerical solution of the Euler equations by ﬁnite volume methods using Runge-Kutta time stepping schemes. AIAA Paper 81-1259, 1981.

[96] P.D. Lax. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. SIAM Regional Series on Appl. Math., 11, 1973.

[85] A. Jameson, W. Schmidt, and E. Turkel. Numerical solutions of the Euler equations by ﬁnite volume methods with Runge-Kutta time stepping schemes. AIAA paper 81-1259, January 1981.

[97] P.D. Lax and B. Wendroﬀ. Systems of conservation laws. Comm. Pure. Appl. Math., 13:217–237, l960.

[86] A. Jameson and E. Turkel. Implicit schemes and LU decompositions. Math. Comp., 37:385–397, 1981.

[98] B. Van Leer. Towards the ultimate conservative diﬀerence scheme. II. Monotonicity and conservation combined in a second order scheme. J. Comp. Phys., 14:361–370, 1974.

[87] M. Jayaram and A. Jameson. Multigrid solution of the Navier-Stokes equations for ﬂow over wings. AIAA paper 88-0705, AIAA 26th Aerospace Sciences Meeting, Reno, Nevada, January 1988.

[99] B. Van Leer. Towards the ultimate conservative diﬀerence scheme. III upstream-centered ﬁnite-diﬀerence schemes for ideal compressible ﬂow. J. Comp. Phys., 23:263–275, 1975. 34

[100] B. Van Leer. Flux vector splitting for the Euler equations. In E. Krause, editor, Proceedings of the 8th International Conference on Numerical Methods in Fluid Dynamics, pages 507–512, Aachen, 1982.

[111] R. L¨ohner and P. Parikh. Generation of threedimensional unstructured grids by the advancing front method. AIAA Paper 88-0515, Reno, NV, January 1988. [112] R.W. MacCormack and A.J. Paullay. Computational eﬃciency achieved by time splitting of ﬁnite diﬀerence operators. AIAA Paper 72154, 1972.

[101] B. Van Leer. Progress in multi-dimensional upwind diﬀerencing. In M. Napolitano and F. Solbetta, editors, Proc. 13th International Conference on Numerical Methods in Fluid Dynamics, pages 1–26, Rome, July 1992. Springer Verlag, 1993.

[113] A. Majda and S. Osher. Numerical viscosity and the entropy condition. Comm. on Pure Appl. Math., 32:797–838, 1979.

[102] B. Van Leer, W. T. Lee, and P. L. Roe. Characteristic time stepping or local preconditioning of the Euler equations. AIAA paper 911552, AIAA 10th Computational Fluid Dynamics Conference, Honolulu, Hawaii, June 1991.

[114] L. Martinelli and A. Jameson. Validation of a multigrid method for the Reynolds averaged equations. AIAA paper 88-0414, 1988. [115] L. Martinelli, A. Jameson, and E. Malfa. Numerical simulation of three-dimensional vortex ﬂows over delta wing conﬁgurations. In M. Napolitano and F. Solbetta, editors, Proc. 13th International Confrence on Numerical Methods in Fluid Dynamics, pages 534–538, Rome, Italy, July 1992. Springer Verlag, 1993.

[103] D. Lefebre, J.Peraire, and K.Morgan. Finite element least squares solutions of the Euler equations using linear and quadratic approximations. Int. J. Comp. Fluid Dynamics, 1:1– 23, 1993.

[116] L. Martinelli and V. Yakhot. RNG-based turbulence transport approximations with applications to transonic ﬂows. AIAA Paper 891950, AIAA 9th Computational Fluid Dynamics Conference, Buﬀalo, NY, June 1989.

[104] S.K. Lele. Compact ﬁnite diﬀerence schemes with spectral-like resolution. CTR Manuscript 107, 1990. [105] J.L. Lions. Optimal Control of Systems Governed by Partial Diﬀerential Equations. Springer-Verlag, New York, 1971. Translated by S.K. Mitter.

[117] D.J. Mavriplis and A. Jameson. Multigrid solution of the Navier-Stokes equations on triangular meshes. AIAA Journal, 28(8):1415– 1425, August 1990.

[106] M-S. Liou and C.J. Steﬀen. A new ﬂux splitting scheme. J. Comp. Phys., 107:23–39, 1993.

[118] D.J. Mavriplis and L. Martinelli. Multigrid solution of compressible turbulent ﬂow on unstructured meshes using a two-equation model. AIAA Paper 91-0237, January 1991.

[107] F. Liu and A. Jameson. Multigrid NavierStokes calculations for three-dimensional cascades. AIAA paper 92-0190, AIAA 30th Aerospace Sciences Meeting, Reno, Nevada, January 1992.

[119] N. D. Melson, M. D. Sanetrik, and H. L. Atkins. Time-accurate Navier-Stokes calculations with multigrid acceleration. In Proceedings of the Sixth Copper Mountain Conference on Multigrid Methods, Copper Mountain, April 1993.

[108] R. L¨ohner and D. Martin. An implicit lineletbased solver for incompressible ﬂows. AIAA paper 92-0668, AIAA 30th Aerospace Sciences Meeting, Reno, NV, January 1992.

[120] J.E. Melton, S.A. Pandya, and J.L. Steger. 3D Euler ﬂow solutions using unstructured Cartesian and prismatic grids. AIAA Paper 93-0331, Reno, NV, January 1993.

[109] R. L¨ohner, K. Morgan, and J. Peraire. Improved adaptive reﬁnement strategies for the ﬁnite element aerodynamic conﬁgurations. AIAA Paper 86-0499, AIAA 24th Aerospace Sciences Meeting, Reno, January 1986.

[121] F. Menter. Zonal two-equation k-ω turbulence models for aerodynamic ﬂows. AIAA Paper 93-2906, AIAA 24th Fluid Dynamics Meeting, Orlando, July 1993.

[110] R. L¨ohner, K. Morgan, J. Peraire, and O.C. Zienkiewicz. Finite element methods for high speed ﬂows. In Proc. AIAA 7th Computational Fluid Dynamics Conference, Cincinnati, OH, 1985. AIAA Paper 85-1531.

[122] E.M. Murman. Analysis of embedded shock waves calculated by relaxation methods. AIAA Journal, 12:626–633, 1974. 35

[123] E.M. Murman and J.D. Cole. Calculation of plane steady transonic ﬂows. AIAA Journal, 9:114–121, 1971.

[136] K. Prendergast and K. Xu. Numerical hydrodynamics from gas kinetic theory. J. Comp. Phys., 109:53–66, November 1993.

[124] R.H. Ni. A multiple grid scheme for solving the Euler equations. AIAA Journal, 20:1565– 1571, 1982.

[137] T.H. Pulliam and J.L. Steger. Implicit ﬁnite diﬀerence simulations of three-dimensional compressible ﬂow. AIAA Journal, 18:159–167, 1980.

[125] S. Obayashi and K. Kuwakara. LU factorization of an implicit scheme for the compressible Navier-Stokes equations. AIAA Paper 841670, AIAA 17th Fluid Dynamics and Plasma Dynamics Conference, Snowmass, June 1984.

[138] J.J. Quirk. An alternative to unstructured grids for computing gas dynamics ﬂows about arbitrarily complex two-dimensional bodies. ICASE Report 92-7, Hampton, VA, February 1992.

[126] J.T. Oden, L. Demkowicz, T. Liszka, and W. Rachowicz. h-p adaptive ﬁnite element methods for compressible and incompressible ﬂows. In S. L. Venneri A. K. Noor, editor, Proceedings of the Symposium on Computational Technology on Flight Vehicles, pages 523–534, Washington, D.C., November 1990. Pergamon.

[139] R. Radespiel, C. Rossow, and R.C. Swanson. An eﬃcient cell-vertex multigrid scheme for the three-dimensional Navier-Stokes equations. In Proc. AIAA 9th Computational Fluid Dynamics Conference, pages 249–260, Buﬀalo, NY, June 1989. AIAA Paper 89-1953-CP. [140] M.M. Rai and P. Moin. Direct numerical simulation of transition and turbulence in a spatially evolving boundary layer. AIAA Paper 91-1607 CP, AIAA 10th Computational Fluid Dynamics Conference, Honolulu, HI, June 1991.

[127] S. Orszag and D. Gottlieb. Numerical analysis of spectral methods. SIAM Regional Series on Appl. Math., 26, 1977. [128] S. Osher. Riemann solvers, the entropy condition, and diﬀerence approximations. SIAM J. Num. Anal., 121:217–235, 1984.

[141] S. V. Rao and S. M. Deshpande. A class of eﬃcient kinetic upwind methods for compressible ﬂows. Report 91 FM 11, Indian Institute of Science, 1991.

[129] S. Osher and S. Chakravarthy. High resolution schemes and the entropy condition. SIAM J. Num Anal., 21:955–984, 1984.

[142] J. Reuther and A. Jameson. Aerodynamic shape optimization of wing and wing-body conﬁgurations using control theory. AIAA paper 95-0213, 33rd Aerospace Sciences Meeting and Exibit, Reno, Nevada, January 1995.

[130] S. Osher and F. Solomon. Upwind diﬀerence schemes for hyperbolic systems of conservation laws. Math. Comp., 38:339–374, 1982. [131] S. Osher and E. Tadmor. On the convergence of diﬀerence approximations to scalar conservation laws. Math. Comp., 50:19–51, 1988.

[143] J. Reuther and A. Jameson. Aerodynamic shape optimization of wing and wing-body conﬁgurations using control theory. AIAA paper 95-0123, AIAA 33rd Aerospace Sciences Meeting, Reno, Nevada, January 1995.

[132] H. Paill`ere and H. Deconinck. A review of multi-dimensional upwind residual distribution schemes for the euler equations. To appear in CFD Review, 1995.

[144] H. Rieger and A. Jameson. Solution of steady three-dimensional compressible Euler and Navier-Stokes equations by and implicit LU scheme. AIAA paper 88-0619, AIAA 26th Aerospace Sciences Meeting, Reno, Nevada, January 1988.

[133] N. A. Pierce and M. B. Giles. Preconditioning on stretched meshes. Report 95/10 1995, Oxford University Computing Laboratory, Oxford, December 1995. [134] O. Pironneau. Optimal Shape Design for Elliptic Systems. Springer-Verlag, New York, 1984.

[145] A. Rizzi and L.E. Eriksson. Computation of ﬂow around wings based on the Euler equations. J. Fluid Mech., 148:45–71, 1984.

[135] K.G. Powell and B. van Leer. A genuinely multidimensional upwind cell-vertex scheme for the Euler equations. AIAA Paper 89-0095, AIAA 27th Aerospace Sciences Meeting, Reno, January 1989.

[146] P.L. Roe. Approximate Riemann solvers, parameter vectors, and diﬀerence schemes. J. Comp. Phys., 43:357–372, 1981. 36

[147] P.E. Rubbert and G.R. Saaris. A general three-dimensional potential ﬂow method applied to V/STOL aerodynamics. SAE Paper 680304, 1968.

[158] R.L. Sorenson. Three-dimensional elliptic grid generation for an F-16. In J.L. Steger and J.F. Thompson, editors, Three-Dimensional Grid Generation for Complex Conﬁgurations: Recent Progress, 1988. AGARDograph.

[148] C.L. Rumsey and V.N. Vatsa. A comparison of the predictive capabilities of several turbulence models using upwind and centered - difference computer codes. AIAA Paper 93-0192, AIAA 31st Aerospace Sciences Meeting, Reno, January 1993.

[159] P. Spalart and S. Allmaras. A one-equation turbulent model for aerodynamic ﬂows. AIAA Paper 92-0439, AIAA 30th Aerospace Sciences Meeting, Reno, NV, January 1992. [160] C.G. Speziale, E.C. Anderson, and R. Abid. A critical evaluation of two-equation models for near wall turbulence. AIAA Paper 90-1481, June 1990. ICASE Report 90-46.

[149] S.S. Samant, J.E. Bussoletti, F.T. Johnson, R.H. Burkhart, B.L. Everson, R.G. Melvin, D.P. Young, L.L. Erickson, and M.D. Madson. TRANAIR: A computer code for transonic analyses of arbitrary conﬁgurations. AIAA Paper 87-0034, 1987.

[161] J.L. Steger and D.S. Chaussee. Generation of body-ﬁtted coordinates using hyperbolic partial diﬀerential equations. SIAM J. Sci. and Stat. Comp., 1:431–437, 1980.

[150] K. Sawada and S. Takanashi. A numerical investigation on wing/nacelle interferences of USB conﬁguration. In Proceedings AIAA 25th Aerospace Sciences Meeting, Reno, NV, 1987. AIAA paper 87-0455.

[162] J.L. Steger and R.F. Warming. Flux vector splitting of the inviscid gas dynamic equations with applications to ﬁnite diﬀerence methods. J. Comp. Phys., 40:263–293, 1981.

[151] C.W. Shu and S. Osher. Eﬃcient implementation of essentially non-oscillatory shockcapturing schemes. J. Comp. Phys., 77:439– 471, 1988.

[163] B. Stouﬄet, J. Periaux, F. Fezoui, and A. Dervieux. Numerical simulation of 3-D hypersonic Euler ﬂows around space vehicles using adapted ﬁnite elements. AIAA paper 870560, January 1987.

[152] C.W. Shu and S. Osher. Eﬃcient implementation of essentially non-oscillatory shockcapturing schemes II. J. Comp. Phys., 83:32– 78, 1989.

[164] R.C. Swanson and E. Turkel. On centraldiﬀerence and upwind schemes. J. Comp. Phys., 101:297–306, 1992.

[153] C.W. Shu, T.A. Zang, G. Erlebacher, D. Whitaker, and S. Osher. High-order ENO schemes applied to two- and three-dimensional compressible ﬂow. Appl. Num. Math., 9:45–71, 1992.

[165] P.K. Sweby. High resolution schemes using ﬂux limiters for hyperbolic conservation laws. SIAM J. Num. Anal., 21:995–1011, 1984. [166] S. Ta’asan, G. Kuruvila, and M. D. Salas. Aerodynamic design and optimization in one shot. AIAA paper 92-005, 30th Aerospace Sciences Meeting and Exibit, Reno, Nevada, January 1992.

[154] B. R. Smith. A near wall model for the k−l two equation turbulence model. AIAA paper 942386, 25th AIAA Fluid Dynamics Conference, Colorado Springs, CO, June 1994. [155] L.M. Smith and W.C. Reynolds. On the Yakhot-Orszag renormalization group for deriving turbulence statistics and models. Phys. Fluids A, 4:364–390, 1992.

[167] Shlomo Ta’asan. Canonical forms of multidimensional, steady inviscid ﬂows. ICASE 93-34, Institute for Computer Applications in Science and Engineering, Hampton, VA, 1993.

[156] R.E. Smith. Three-dimensional algebraic mesh generation. In Proc. AIAA 6th Computational Fluid Dynamics Conference, Danvers, MA, 1983. AIAA Paper 83-1904.

[168] E. Tadmor. Numerical viscosity and the entropy condition for conservative diﬀerence schemes. Math. Comp., 32:369–382, 1984.

[157] R.L. Sorenson. Elliptic generation of compressible three-dimensional grids about realistic aircraft. In J. Hauser and C. Taylor, editors, International Conference on Numerical Grid Generation in Computational Fluid Dynamics, Landshut, F. R. G., 1986.

[169] J.F. Thompson, F.C. Thames, and C.W. Mastin. Automatic numerical generation of body-ﬁtted curvilinear coordinate system for ﬁeld containing any number of arbitrary twodimensional bodies. J. Comp. Phys., 15:299– 319, 1974. 37

[170] J.F. Thompson, Z.U.A. Warsi, and C.W. Mastin. Boundary-ﬁtted coordinate systems for numerical solution of partial diﬀerential equations: A review. J. Comp. Phys., 47:1– 108, 1982.

[182] H.C. Yee. On symmetric and upwind TVD schemes. In Proc. 6th GAMM Conference on Numerical Methods in Fluid Mechanics, Gottingen, September 1985. [183] S. Yoon and A. Jameson. Lower-upper Symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations. AIAA Paper 870600, AIAA 25th Aerospace Sciences Meeting, Reno, January 1987.

[171] E. Turkel. Preconditioned methods for solving the incompressible and low speed equations. J. Comp. Phys., 72:277–298, 1987. [172] V. Venkatakrishnan. Newton solution of inviscid and viscous problems. AIAA paper 880413, January 1988. [173] V. Venkatakrishnan. Convergence to steady state solutions of the Euler equations on unstructured grids with limiters. AIAA paper 930880, AIAA 31st Aerospace Sciences Meeting, Reno, Nevada, January 1993. [174] G. Vorono¨ı. Nouvelles applications des parametres continus `a la th´eorie des formes quadratiques. Deuxi`eme m´emoire: Recherches sur les parallelloedres primitifs. J. Reine Angew. Math., 134:198–287, 1908. [175] Y. Wada and M-S. Liou. A ﬂux splitting scheme with high-resolution and robustness for discontinuities. AIAA paper 940083, AIAA 32nd Aerospace Sciences Meeting, Reno, Nevada, January 1994. [176] N.P. Weatherill and C.A. Forsey. Grid generation and ﬂow calculations for aircraft geometries. J. Aircraft, 22:855–860, 1985. [177] D.C. Wilcox. A half a century historical review of the k-ω model. AIAA Paper 91-0615, AIAA 29th Aerospace Sciences Meeting, Reno, NV, January 1991. [178] F. Woodward. An improved method for the aerodynamic analysis of wing-body-tail conﬁgurations in subsonic and supersonic ﬂow, Part 1 — theory and application. NASA CR 2228 Pt. 1, May 1973. [179] P. Woodward and P. Colella. The numerical simulation of two-dimensional ﬂuid ﬂow with strong shocks. J. Comp. Phys., 54:115–173, 1984. [180] K. Xu, L. Martinelli, and A. Jameson. Gaskinetic ﬁnite volume methods, ﬂux-vector splitting and artiﬁcial diﬀusion. J. Comp. Phys., 120:48–65, 1995. [181] V. Yakhot and S.A. Orszag. Renormalization group analysis of turbulence. I. Basic theory. J. Sci. Comp., 1:3–51, 1986. 38

Figure 18: Navier-Stokes Predictions for the F-18 Wing-Fuselage at Large Incidence

39

1.40 1.00

2.00 0.00

0.80

-2.00

0.60

-4.00

0.40

-6.00

+

0.20

+

++++++++++ ++++ ++++ +++ +++ +++ +++ ++ + + + + +

-8.00

+

+

0.40

Log(Error)

+++++++ + + + + + ++++ ++++++++ ++++++ +++++ ++++ ++++ +++ +++ +++ +++ +++ +++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + ++ + + ++ ++ ++ ++ + + ++ ++ ++ ++ + ++ ++ ++ ++ ++ + + ++ +++ +++ + + ++++

+

1.20

4.00

-2.00 -1.60 -1.20 -0.80 -0.40 0.00

+

0.00

-10.00

+ + + + + + +

-0.20

-12.00

0.80

+

1.20

Cp

++++++ +++++ ++++ ++++ ++++ + + + + ++ +++++++ +++++++++ +++++++++++++++ +++++++++++++++++++++++++++++++++ ++ + + + + + + + + + + + + + + + + + + +

0.00

50.00

100.00

150.00

200.00

Work

19a: Cp after 25 Cycles. Cl = 1.1312, Cd = 0.0469.

19b: Convergence.

Figure 19: RAE-2822 Airfoil at Mach 0.750 and α = 3.0◦ H-CUSP Scheme.

40

250.00

300.00

1.40

0.00

1.00

2.00

1.20

4.00

-2.00 -1.60

0.80

-2.00

0.60

-4.00

-0.20

-12.00

0.00

-10.00

0.20

-8.00

0.40

-6.00

Log(Error)

+++++++++ +++++++ +++++ +++++ ++++ ++++ ++++ +++ +++ +++ +++ +++ +++ ++ + ++ ++ ++ + ++ ++ ++ ++ ++ + ++ ++ ++ ++ ++ + ++ ++ ++ ++ ++ + + ++ ++ ++ ++ +++ + + + +++ +++ ++++ + + + + ++++

+++++ + + + + + + + + + + + + + + + + + + + ++

-1.20 -0.80 -0.40

Cp

0.00 0.40 0.80 1.20

+++++++++++++ ++++++++++ +++++++ ++++++ +++++ +++++ + + + ++++ ++++ ++++ +++ +++ + + + +++ +++ ++ ++ ++ + + + + + + + + + + + + + ++++++++++++++++ +++++ + ++++ + +++ + +++ + +++ ++ + ++ + ++ + + + + + + + + + + + + + + + + + + +

0.00

50.00

100.00

150.00

200.00

250.00

300.00

Work

20a: Cp after 35 Cycles. Cl = 0.3654, Cd = 0.0232.

20b: Convergence.

1.40 1.20 1.00

0.00

0.80 0.60

-4.00

-0.20

-12.00

0.00

-10.00

0.20

-8.00

0.40

-6.00

Log(Error)

+++++++++ +++++++ +++++ +++++ ++++ ++++ ++++ +++ +++ +++ +++ +++ +++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + +++ ++ + +++ ++ ++ ++ ++ ++ + ++ +++ +++ +++ +++ + + + + + + ++++++

-2.00

+++++++++++++++++ ++++++++++++++ +++++++++ +++++++ ++++++ +++++ + + + + ++ ++++ ++++ ++++ ++++ +++ + + +++ +++ +++ +++ ++ + + ++ + + + + + + + + + + + + + ++++ + +++ +++ + +++ + +++ + ++ + ++ + + + + + + + + + + + + + + + + + + + +

++++ + + + + + + + + + + + + + + + + + + + + ++

-0.40 0.00 0.40 0.80 1.20

Cp

-0.80

-1.20

2.00

-1.60

4.00

-2.00

Figure 20: NACA-0012 Airfoil at Mach 0.800 and α = 1.25◦ H-CUSP Scheme.

0.00

50.00

100.00

150.00

200.00

Work

21a: Cp after 35 Cycles. Cl = 0.3861, Cd = 0.0582.

21b: Convergence.

Figure 21: NACA-0012 Airfoil at Mach 0.850 and α = 1.0◦ H-CUSP Scheme. 41

250.00

300.00

-2.00 -1.20

-1.60

-2.00 -1.60 -1.20

+++ + ++ + +

+++

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+ +

1.20

0.80

+ + +

++ + +

++ + +

+

1.20

-2.00 -0.80

+++

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

++

+ +

++

+

+++ + + + + + +

+

+

+

+

+

+

+ + ++ +++++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + + + + +

0.00

+

+

+

+

+

+

+

+

+

+

+

++

+

++

0.40

+ ++ +++++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + +

++++ + +

+

++ + ++ ++ + +++ ++

+

+

+

++++ ++ + + + +++++++ ++++ + + +

+ -0.40

+

+ +++ +++ +++ +++++ ++ ++

Cp

++++ ++ + + ++ + ++ + + + + +

-1.20

-1.60

-2.00

22b: 31.25% Span. Cl = 0.3139, Cd = 0.0159.

+

+

+

+

+

0.80

++ + +

+

+ +

1.20

1.20

+ +

++ +

0.80

+

+ 0.40

+

+

+

-1.60 -1.20 -0.80 -0.40

Cp

++

+

+

+

22a: 12.50% Span. Cl = 0.2933, Cd = 0.0274.

0.00

+++++

+

+

+

+

+

+

++ +

+ ++ +++++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + + + + +

+

+++

+

+

+

0.40

-0.40

Cp +

+

0.00

++++

+

++ + +++ +++ ++++++++++++ ++ ++

+

+

+

0.40

+

+

+

0.80

+ +

+

+ + ++ + + +

+ + + + +++++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + +

-0.40

Cp

+ ++ + +++ + +++ +++++ ++++++++++++++++ +

+

-0.80

+

+

0.00

-0.80

+

+

22c: 50.00% Span. Cl = 0.3262, Cd = 0.0089.

22d: 68.75% Span. Cl = 0.3195, Cd = 0.0026.

Figure 22: Onera M6 Wing. Mach 0.840, Angle of Attack 3.06◦ , 192×32×48 Mesh. CL = 0.3041, CD = 0.0131. H-CUSP scheme. 42

23a: Initial Wing CL = 0.5001, CD = 0.0210, α = −1.672◦

23b: 40 Design Iterations CL = 0.5000, CD = 0.0112, α = −0.283◦

Figure 23: Swept Wing Design Case (1), M = 0.85, Fixed Lift Mode.Drag Reduction at CL = .5.

43

UPPER SURFACE PRESSURE

UPPER SURFACE PRESSURE

24a: Initial Wing CL = 0.5001, CD = 0.0210, α = −1.672◦

24b: 40 Design Iterations CL = 0.5000, CD = 0.0112, α = −0.283◦

Figure 24: Swept Wing Design Case (1), M = 0.85, Fixed Lift Mode.Drag Reduction at CL = 0.50

44

-2.00 -1.20

-1.60

-2.00 -1.60 -1.20

+

-.80

+

Cp

.00

+

++ ++++++++++++ +++ +++ ++ +++++++++ ++ + + ++ + + +

++

+

++

+

+

+

++

+

+

+

+

+

+

++

+

+ + + + + + + + +

++ +

+

+ + +

+ + + +

+

+

+

.40

+

+

+

+

+ +

+

+

++

+

+ + ++ ++

+ ++ + +++ ++++++++++++++++ ++ + + + ++ ++ + ++ + + + + + +

+ ++ +

++

+ + ++ ++++ +++ +++ ++ ++ ++ ++ + ++ + + + + + + ++ + + + + + + + + + + + + +

+

+

+++ +++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + +

+

+

+++

-.40

-.80 -.40

Cp

.00 .40

++

+ ++

+

.80

.80

+

+

+

+

1.20

1.20

++ +

+ + +

25a: span station z = 0.00

-2.00 -1.60 -1.20 -.80 -.40

+ + + ++ ++++ +++ +++ ++ ++ ++ ++ + ++ +++ + + ++ + + + + + + + + + + + + + + + +

.00

+ +

+

+ + + + +

+

+

+

+

+ + +

.80 1.20

1.20

.80

+ +

+ +

+ + +

25c: span station z = 0.625

25d: span station z = 0.937

Figure 25: FLO67 solution for initial wing.M = 0.85, CL = 0.4997, CD = 0.0207, α = −1.970◦. 45

+

+

+

.40

+

+ + + + + + +

+

+

+

+

+

+

+

+ + + + + + +

+

.40

+

+

+

+

+

++ ++ + +++ + ++ ++ ++ ++ +++ +++ +++++++++++++++ ++++ + + + + + ++ + +

+ +

+

.00

Cp

++ +++ ++++ ++ + + + + +++ +++ + ++ + ++ +++ + + + ++ +++++++++ ++ + + + + ++++

+ + + ++ ++++ +++ +++ ++ ++ ++ ++ + ++ ++ +++ ++ + + + + + + + + + + + + + + + + +

-.40

Cp

-.80

-1.20

-1.60

-2.00

25b: span station z = 0.312

-2.00

++

++

+++

+ +

+ +

+ +

+

+ +

+ + + +

+ + +

+

+

++ + + + + + + + ++

.40

+

+

+

+

+

+

+++ +++ ++ ++ ++ ++ ++ ++ + ++ +++ + + + + + + + + + + + + + + + +

+++ +

+

-1.60

Cp

.00

+

+

+

+

+

+ +

+++ ++ ++ ++ ++ ++ ++ ++ ++ +++ + + + + + + + + + + + + + + + + + + + +

+ + ++ + ++

+

.40

-.40

-.80

-1.20

-2.00 -1.60 -1.20 -.80 -.40

Cp

.00

++ ++ +++ + + + + + + + + + +++++ + +++++++++++ + ++++++ +++++ + +++ + + + + + + + + +

+ +++++++++ + +++++++++++ ++++ +++ + ++++ ++++ + + +

+

+

+ +

.80

++ + + + + +

.80

+

+ + +

1.20

26a: span station z = 0.00

-2.00 -1.60 -1.20

+

+

+

+

+

.00

+

+

+

+

.40

+

.40

++ + + + + + + ++

+

+

+

+ +

++

++ + + +

++++

+++

+

+

+

+ + +

+ + + + + + + +

+

+

++

+

++

+

++

-.80

++ ++ + +

++++ ++ +++++++++++++++ ++++++ + +++ ++ + ++ + +

-.40

+

Cp

++++++++ + +++++++++ ++++ ++++ + ++++ ++ + +++++ + +

++ + + + + + + ++ +++ +++ ++ ++ ++ ++ ++ ++ + ++ +++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + +

-.40 .00

Cp

-.80

-1.20

-1.60

-2.00

26b: span station z = 0.312

+++ +++ ++ ++ ++ ++ ++ ++ ++ ++ +++ + + + + + + + + + + + + + + + + +

1.20

+

+

+ + +

.80 1.20

1.20

.80

+ +

26c: span station z = 0.625

26d: span station z = 0.937

Figure 26: FLO67 check on redesigned wing.M = 0.85, CL = 0.4992, CD = 0.0094, α = −0.300◦. 46

-2.00

+

+

+

+ +

+

+ +

+

+

+ + +

+

+ +

+ + + +

+

.40

+

+

+

+

+

+

+

+++

+

+

+

+

+

+ +

.80

+

++ + + + + +

.80

+

+ + + + + + + + ++ +

.00

+ +

+

+ + +++ +++ ++ ++ ++ ++ ++ ++ + ++ +++ + + + + + + + + + + + + + + +

-1.20

+

+++ ++ ++ ++ ++ ++ ++ ++ ++ +++ + + + + + + + + + + + + + + + + + + + +

+ + ++ + ++

Cp

++++ + + + + + + + +

-.40

-.80

-1.60

-2.00 -1.60 -1.20 -.80 -.40

Cp

.00 .40

+++ ++ +++ ++++++++ ++++++++ ++++++ ++++ + + + + + + + + +

+ ++++++++++++++++ +++++ + +++++ ++++ ++++++ +++ + ++++ + +

+ + +

1.20

27a: span station z = 0.00

-2.00 -1.60 -1.20

+

+

+

.00

.00

+

+

+

.40

+

.40

++ + + + + + + ++

+

+

+

+ +++ ++ + + +

++++

+++

+

+

+

+ + +

+

+

+

+ + + + + + +

+

+

+

+

++

+

+ ++ + + +

++ + + + + + + ++ +++ +++ ++ ++ ++ ++ ++ ++ + ++ +++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + +

++

++++ ++ +++++++++++++++++++++++++++++ + + +

+

+

-.40

+ +

-.80

+++++++++ + ++++++++++++++ +++++++++++++ ++ + +

Cp

-.40

Cp

-.80

-1.20

-1.60

-2.00

27b: span station z = 0.312

+++ +++ ++ ++ ++ ++ ++ ++ ++ ++ +++ + + + + + + + + + + + + + + + +

1.20

+

+ +

.80 1.20

1.20

.80

+ +

27c: span station z = 0.625

27d: span station z = 0.937

Figure 27: FLO67 check on redesigned wing at a higher Mach number. M = 0.86, CL = 0.4988, CD = 0.0097, α = −0.440◦. 47