Introduction

Panel methods are numerical schemes for solving (the Prandtl-Glauert equation) for linear, inviscid, irrotational flow at subsonic or supersonic free-stream Mach numbers. Currently, panel-method codes are the only codes commonly in use that are sufficiently developed for routinely analyzing the complex geometries of realistic aircraft. The objective of this report is to give the reader some idea of what panel methods can and cannot do, to describe their common roots, to describe the difference between various specific implementations, and to show some example applications. In addition, recent progress in solving nonlinear transonic flow problems by combining portions of panel-method technology with other numerical techniques, is described. This material is followed by three appendixes that contain additional details: appendix A describes how the basic ideas common to all panel methods are actually implemented in a specific code; appendix B shows how to evaluate some of the integrals that arise in the influence-coefficient equations; and appendix C discusses the so-called finite part of improper integrals.

There are fundamental analytic solutions to the Prandtl-Glauert equation known as source, doublet, and vorticity singularities. Panel methods are based on the principle of superimposing surface distributions of these singularities over small quadrilateral portions, called panels, of the aircraft surface, or to some approximation to the aircraft surface. The resulting distribution of superimposed singularities automaticlaly satisfies the Prandtl-Glauert equation. To make the solution correspond to the desired geometry, boundary conditions are imposed at discrete points of the panels. (Mathematicians refer to these discrete points as collocation points; panel-method users refer to them as control points.)

Panel codes are often described as being lower-order or higher-order. The term lower-order refers to the use of constant-strength singularity distributions over each panel, and the panels are usually flat. Higher-order codes use something of higher order than constant, for example, a linear or quadratic singularity distribution, and sometimes curved panels.

Panel methods were initially developed as lower-order methods for incompressible and subsonic flows (e.g., refs. 1,2; see ref. 3 for a review of panel methods existing through about 1976). The first successful panel method for supersonic flow became available in the mid-1960s (refs. 4,5). This was also a lower-order method, and is variously referred to as the constant-pressure panel method, or the Woodward-Carmichael method.

The panel methods for three-dimensional subsonic flow allowed the actual vehicle surface to be paneled, whereas the Woodward-Carmichael method was more severely restricted in the placement of the panels. For example, the wing was a planar array of panels, the body (fuselage) volume was modeled with a line distribution of source and doublet singularities (resulting in a body-of-revolution) and the body boundary conditions were imposed with a cylindrical "interference" shell of wing-type panels. This representation was later extended to include multiple wing-body components (ref. 6), but was still restricted to the planar panel representation. These two extremes of actual-surface models and mean-surface models (Woodward-Carmichael) are illustrated in figure 1.

The mean-surface model used in the Woodward-Carmichael panel method was a consequence of numerical stability problems that arose in supersonic flow. The constant-strength, elementary horseshoe vortex singularity distribution (producing a constant pressure over each panel) often produced unstable numerical behavior (the solutions would "blow-up") when a panel was inclined to a supersonic flow. The method worked only when all the panels were kept parallel to the free-stream flow. This required that angle of attack, wing thickness, camber, and twist be simulated through the boundary conditions; that is, it was necessary to have the panels generate flow perpendicular to themselves and thereby turn the flow through the desired angles, as is done in classic thin-wing theory. As a consequence of this restricted geometric model, several new approaches to the supersonic problem were pursued in the 1970s.

The first of these was also due to Woodward; it evolved into the series of computer programs known as USSAERO (ref. 7). For fuselage panels, USSAERO uses constant-strength source singularities whose strength distribution varies linearly in the chordwise direction and is constant in the spanwise direction. Although this representation gave an improved modeling capability, numerical problems would still often occur when the wing panels were inclined to supersonic free-stream flow.

Another approach, developed by Morino and his associates, uses a superposition of constant-strength sources and doublets on hyperboloidal panels. The constant-strength doublets produce a velocity field that is identical to that produced by line-vortex elements, having the same strength as the doublet panel, placed head to tail around the panel perimeter (a so-called ring-vortex panel). This method is available in the computer program called SOUSSA; it too is unable to handle the steady supersonic case (ref. 8, pp. 2-20; private communication, L. Morino, Feb. 1981).

The key to eliminating the numerical stability problems associated with supersonic flow, was to use doublet distributions that were continuous over the entire surface of the aircraft. This approach, using quadratic doublet distributions (equivalent to linear vorticity distributions) was first used in the PAN AIR code (refs. 9-14) and its pilot code predecessor (ref. 15). It has since been implemented in the European version of PAN AIR, called HISSS (ref. 16). The continuous-doublet distribution eliminates the appearance of spurious line-vortex terms at the panel edges, which was the cause of the numerical stability problems in the earlier approaches.

Within the limitations of the Prandtl-Glauert equation, the higher-order singularity distributions used in the PAN AIR and HISSS codes allowed the actual-surface paneling models, long in use for subsonic flow, to also be used for supersonic flow. It also had a very beneficial side effect: the numerical solutions turned out to be much less sensitive to the size, shape, and arrangement of the paneling than in earlier methods, including the subsonic-only methods. Partly for this reason, continuous quadratic doublets were incorporated into the subsonic-only MCAERO code (ref. 17). These advantages did not come free however. The higher-order distributions require much more analytic work to derive the influence-coefficient equations, and demand many more arithmetic operations than the simpler lower-order (constant-strength) methods, which results in significantly higher run costs.

It was subsequently discovered that for subsonic flow, setting the pertubation potential to zero at the interior side of panels, in conjunction with the original lower-order singularity distributions, also reduced the solution sensitivity to variations in panel layout. This led to a renewed interest in the lower-order methods, resulting in the VSAERO (refs. 18, 19) and QUADPAN (refs. 20, 21) codes. QUADPAN was later revised to handle the supersonic case by changing its constant-strength doublets to continuous linear doublets.

What Panel Methods Can and Cannot Do

Panel-method-based computer programs are currently the workhorse codes for predicting the aerodynamics of complete configurations. Representative aircraft examples that have been analyzed with panel method codes are shown in figure 2. Although such codes are routinely used to analyze very complicated geometries, they do so at the expense of ignoring much fluid physics. The equation that panel codes solve is the Prandtl-Glauert equation. For steady subsonic flow this equation is usually written as
${\tilde{\nabla}}^{2} ϕ = (1 - M_{\infty}^{2}) ϕ_{x x} + ϕ_{y y} + ϕ_{z z} = 0$
and for supersonic flow it is sometimes multiplied by -1,
$- {\tilde{\nabla}}^{2} ϕ = (1 - M_{\infty}^{2}) ϕ_{x x} + ϕ_{y y} + ϕ_{z z} = 0$
where $M_{\infty}$ is the free-stream Mach number and $ϕ$ is the pertubation velocity potential.

For subsonic flee-stream flow, equations (1) and (2) are elliptic being similar to Laplace's equation. Such equation types have the property that any disturbance at some point is felt everywhere in the flow field (although the effect usually dies out rapidly with distance). For supersonic free-stream flow the equations are hyperbolic, with the x-derivative term behaving like time in the wave equation. Solutions for the supersonic case are fundamentally different, disturbances having restricted zones of influence (or in Von Karman's words, zones of "silence", or "forbidden signals"; ref. 22). The disturbances only propogate downstream, along rays defined by the Mach cones (characteristic surfaces), reflecting off downstream geometry and interfering in a wave-like manner with other disturbances.

The Prandtl-Glauert equation is the simplest form of the fluid-flow equations that contain compressibility effects (i.e., the effect of Mach number on fluid density). It is obtained from the more general Navier-Stokes equation by (1) neglecting all the viscous and heat-transfer terms; (2) assuming that the flow is irrotational, thereby admitting the introduction of a velocity potential; and (3) discarding all nonlinear terms. This restricts the flow to be inviscid, irrotational, and linear. Often, the flow is also assumed to be steady. Physically, these restrictions mean that important flow behavior such as separation, skin-friction drag, and transonic shocks are not predicted with panel methods. Items that are predicted include drag-due-to-lift (often called induced drag for subsonic flow, and vortex drag for supersonic flow), and wave drag.

Wave drag is predicted because the Prandtl-Glauert equation admits solutions that approximate the weak-shock solutions of shock-expansion theory (ref. 23, pp. 215, 216). A simple example is the supersonic flow over a thin wedge (fig. 3(b)). For small wedge (deflection) angles, the shock is attached at the wedge leading edge, forms at an angle very nearly to that of the Mach angle, and the flow remains supersonic on the downstream side of the shock. The limiting case for these weak shocks, in which the shocks form at exactly the Mach angle, is pre-dicted by the Prandtl-Glauert equation.

The absence of any explicit viscous effects causes subsonic flow solutions to be non-unique unless a Kutta condition at sharp trailing edges is somehow imposed (ref. 24, pp. 80, 81). This is done with the addition of some type of wake panels that trail downstream from lifting-surface trailing edges (fig. 3(a)), causing the flow to separate smoothly from these edges and allowing the potential to jump (be discontinuous) across the wake. Most panel methods require the user to assume the shape and position of the wakes. For a simple wing body this poses no difficulty, the wake position being relatively unimportant. However, for multiple-lifting-surface config-urations, the wake placement is important since it affects the flow experienced by downstream geometry. A few codes iteratively solve for the wake shape and location.

Because panel methods are able to treat complete configurations, they have often been used in combination with other methods to approximately account for additional physics neglected by the Prandtl-Glauert equation. One fairly common practice is to include the presence of the wing boundary layer (ref. 25). The basic idea is to use the pressure distribution from the panel-code solution as input to a boundary-layer code and compute the displace-ment thickness. This incremental thickness is then repre-sented in a second run of the panel code. This is usually done in one of two ways, as illustrated in figure 4 (ref. 25). The first is to actually recompute the wing surface coordinates and the new wing-body intersection by adding the displacement thickness to the actual wing geometry. An alternative approach is to use "blowing," in which the source strengths of the wing panels are adjusted such that each panel ejects (or sucks) enough fluid to cause the resultant flow field to be approximately displaced by the displacement thickness. For either approach, the resulting change in actual or apparent wing shape has two effects: it reduces the effective camber of a cambered wing and it increases the wing thickness. For a specified angle of attack, the primary aerodynamic effect of these changes is a reduced lift owing to the reduced camber. The second, but usually less important effect, is a slight increase in lift owing to the increased wing thickness.

Another example is the coupling of panel codes to propulsion codes. In reference 26, the PAN AIR code is coupled to a parabolized Navier-Stokes propulsion code. The purpose was to account for the viscous, high-energy, exhaust-flow effect on the aerodynamic flow about the complete aircraft.

Panel-method codes have also been built to model the flow separation that occurs off highly swept wings with sharp leading edges (refs. 27, 28). In these codes, wake panels emanate from the wing leading edge, as well as from the trailing edge (fig. 5). Iterative techniques are used to solve for the correct shape and position of the leading-edge wake panels. The criteria to be satisfied are (1) that the Kutta condition be enforced and (2) that the entire wake surface be a stream surface (i.e., no flow crosses it, and it supports no pressure jump).

Common Roots of Panel Methods

As indicated in section 2, panel methods rely on sur-face distributions of sources, doublets, and vorticity. We will see later that doublet and vortex distributions are related. Since surface vorticity is a vector and a doublet is a scalar, it is often easier to work with doublets than with vorticity, and then compute the vorticity from the doublet-strength distribution. Most higher-order panel-method codes take this approach.

It can be verified by direct substitution, that the following expressions, called unit point sources and doublets, respectively, satisfy the Prandtl-Glauert equation (eqs. (1) or (2)). $ϕ_{P}^{S} ({\vec{X}}_{Q}) = \frac{- 1}{R ({\vec{X}}_{P}, {\vec{X}}_{Q})}$

In these expressions, point P is the influenced point in space having coordinates _p = (xp, yp, zp), and point Q is the influencing point gQ = (xQ, yQ,ZQ) at which the unit point source or doublet is located (see fig. 6). There is an elemental area dSQ associated with the doublet, and the doublet axis is normal to this area. (Recall from elemen-tary fluid mechanics that a doublet can be thought of as a source-sink pair approaching each other along an axis. This definition of a doublet produces the same result as eq. (4).) The subscript Q on the scaled gradient operator means that the derivatives are to be taken with respect to the coordinates of point Q, not point P.

For incompressible flow, R becomes simply the geo-metric distance between the two points P and Q. Equation (3) then tells us that a point source at Q produces a distur-bance at P that diminishes inversely as the distance between the two points. The meaning of equation (4) is not so obvious until one works out the expression indi-cated by the dot product. If one chooses the xyz coordi-nate system at point Q as shown in figure 6, then the unit normal fi equals the unit vector k and equation (4) becomes simply This form clearly shows the directional properties of a point doublet and reveals that a doublet disturbance dies off at least as rapidly as the inverse of the distance squared.

Since the Prandtl-Glauert equation is a linear partial differential equation, sums of the source and doublet solutions are also solutions. Thus, panel methods are usu-ally thought of as superposition methods, and, hence, are restricted to linear problems. There is a more general approach, however, that, while containing superposition as a special ease, can also be used to solve nonlinear prob-lems. In section 6 we will take a look at how panel-method technology can be combined with other tech-niques to solve the nonlinear full-potential equation, so it is advantageous to look at this more general approach, known as Green's third theorem (ref. 29, p. 21, eq. (7)).

In reference 29, the derivation corresponds to incompressible potential flow; in reference 30, this is generalized to the compressible case. The result is the following identity: (eq7) where (eq8).

In the above equations (sigma) is the source strength and (mu) is the doublet strength at any point Q, on the surface S, which in our case will be all (for subsonic flow) or part (for supersonic flow) of the aircraft surface and wakes. These strengths are equal to jumps (discontinuities) across the panels of certain flow properties. The source strength equals the jump in the normal component of the mass-flux vector ft. The doublet strength equals the jump in potential and the gradient of the doublet strength equals the jump in tangential component of velocity. he values of these strengths are the (as yet unknown) amplitudes of the source and doublet singularity solutions appearing in equations (3) and (4). Here, these source and doublet solutions, when multiplied by a constant k-1, are denoted as Ka and K_t, respectively (K is used to denote that the singularities are called kernels). For M** < 1, k = 4r¢, and S is the entire surface of the aircraft and wake(s). For M** > 1, k = 2_t, and S is that portion of the aircraft surface and wake(s) that lies in the upstream Mach cone emanating from the influenced point P.

Equation (7) shows us that the velocity potential at point P is related to the source and doublet distributions on S, and to the spatial distribution of V2 0 in the volume V bounded (wetted by) both sides of S. If d_p is con-structed according to the surface integral terms in equa-tion (7), that is, (eq 9) then, because equation (7) is an identity, it follows that (eq 10). Since K_{\sigma} is a function of the arbitrary point P, (\tildegrad^{2}\phi) must be zero. Thus, construction of \phi_{P} according to equation (9) implies that equations (1) and (2), the Prandtl-Glauert equation, has been satisfied throughout V.

Equation (9) is the basic starting forumla for panel methods using sources and doublets. If the source and doublet strength distribution is known (we will see how this is done later), then the velocity at point P is obtained from equation (9) by differentiating with respect to the coordinates at P, that is (eq 11). Equations (9) and (11) are used to generate influence-coefficient equations that relate source and doublet strengths at particular points Q on the surface S to the potential and velocity at field points P. The basic idea is to break S into a collection of panels \Sigma and to assume a functional form for \sigma and \mu over each panel. For example, a constant-strength source-doublet panel with index j is given simply by (eq 12a) (eq 12b) where the unknown constants \lambda^{S}_{j} and \lambda^{D}_{j} are called source and doublet singularity parameters, respectively, for panel j, and (\zeta, \eta) are local surface coordinates associated with the panel.

Once the functional form for aj(_,rl) and Bj(_,_) are specified, equations (9) and (11) can be integrated over each panel (a nontrivial task) so that 0p and Wp are expressions involving only the unknown singularity parameters. If P is made a control point (a panel point at which a boundary condition will be imposed) with index i, equations (9) and (11) give the potential and velocity at that point in terms of (as influenced by) the source and doublet disuibutions of the single panel j (see fig. 7). Note that the fixed point P and the variable point Q of the ana-lytic formulation correspond to control point i and panel j, respectively, in the discretized implementation. Summing the effects from all the panels on the aircraft surface gives the potential and velocity at control point i in terms of the total number (N) of singularity parameters.

If \arrow{v}_{ij} denotes the velocity at control point i, owing to the source-doublet distributions at panel j, then the velocity at point i owing to all N panels is (eq 13). If the panel associated with control point i is to be a solid (impermeable) panel represented by a zero normal component of the total velocity, then the boundary condition is (eq 14). Thus, for control point i, we have (eq 15) and since the \arrow{v}_{ij} are known (from eq. (11)) in terms of the N singularity parameters, equation (15) can ultimately be expressed as the single equation (eq 16) where the VIC_{ij} are called velocity influence coefficients and b_{i}=-\arrow{V}_{\infty}\dot\hat{n}_{i}.

Repeating the process for a total of N control points and applying boundary conditions at each of these points leads to the equation (eq 17). In this equation, [AIC] is called the matrix of aerodynamic influence-coefficients, {\lambda} is the vector of unknown singularity parameters, and the elements of {b}, the so-called right-hand-side quantities, are known from the boundary conditions. Solving for the singularity parameters then makes it possible to compute (from eqs. (9) and (11)) the potential and velocity distributions, and hence the pressures acting at the paneled approximation to the surface S. Integrating the pressures and their moments yields the resultant aerodynamic force and moment.

To summarize the basic ideas, a panel method uses the fact that sources and doublets are solutions to the Prandtl-Glauert equation. The numerical procedure is as follows: 1. Break the aircraft surface into an assemblage of panels 2. Use the panels to create source-doublet distributions in terms of singularity parameters {\lambda}, whose values are to be determined 3. Form influence-coefficient expressions for the potential or velocity or both at each of N panel control points owing to the source-doublet distribution of any panel; for supersonic flow the zones of silence must be accounted for 4. Use the influence-coefficient expressions to enforce boundary conditions at N points, giving the NxN system of equations (eq systems) 5. Solve for the singularity parameters {\lambda}, then compute the potential and velocity at any point of interest 6. Knowing the velocity, compute the panel pressure distributions 7. Integrate the pressure distributions to obtain forces and moments The above process, in one form or another, is common to all panel methods. Selecting a way to break the surface S into panels \Sigma, choosing and implementing funtional forms for \sigma and \mu, and evaluating the influence-coefficient surface integrals are the fundamental tasks.

In section 4 we discuss how various codes differ in the selection of I;, and of o(_,T1) and It(_,'q). Section 5 describes various means of specifying boundary condi-tions and other aspects of modeling physical problems. Appendix A gives additional details about how the basic approach described above is actually implemented into the PAN AIR code. This will lead to specific surface inte-grals that must be evaluated to obtain the influence-coefficients. The analytic evaluation of some of these integrals for subsonic flow is given in appendix B.

The reader is warned that the above description has ignored some subtle points that will be addressed later. Specifically, a panel with both source and doublet distri-butions requires two boundary conditions per panel. These boundary conditions generally require some statement about the flow on each side of the panel. (There is an exception in supersonic flow.) The boundary conditions must also produce a well-posed mathematical problem, that is, one with a unique solution; otherwise the AIC matrix will be singular. Also, some of the analytic expres-sions for the potential and velocity can become singular (blow up) when the influenced point is a point of the influencing panel (when i = j); hence use of the name source-doublet "singularities." Boundary conditions will be discussed in section 5, and improper integrals (those having a nonintegrable singularity) will be discussed in appendix C.

DISTINCTIONS BETWEEN PANEL CODES

Having outlined the basic approach used by panel methods, let us now look at various implementations. Dif-ferent ways of approximating the actual aircraft surface with surface panels will be described first, followed by the selection of singularity distributions and their functional forms. Finally, a summary of several specific panel codes will be presented.

Surface Geometry

Current actual-surface panel codes generally allow the user to break the configuration into logical pieces such as forebody, canopy, and the wing upper and lower sur-faces. In the PAN AIR code, such logical pieces are called networks of panels (fig. 8). A component such as a wing can also be broken into several networks. For example, two different models for the F-16XL wing upper surface are shown in figure 9. The 3-network model is appropriate if control-surface deflections are not of interest; the 10-network model allows leading and trailing-edge control surface deflections to be made.

Panel codes typically use arrays of MxN surface grid points to define the comer points of panels, as illus-trated in figure 10. Since the four panel comer points are not, in general, coplanar, a single flat panel cannot be used to connect the four corner points.

The lower-order (constant-strength) codes usually use some kind of flat "average" panel defined by the four comer points. Using flat quadrilateral panels to represent a curved surface causes gaps to exist between the edges of adjacent panels (fig. 11, courtesy of John Hess). This approach begs the natural question about "leaks" through the gaps. Actually, for constant-strength panels, the panels themselves leak rather badly everywhere except at the control points x, where the discrete boundary conditions are imposed. The constant-strength, singularity influence-coefficient-computed velocities at any points on the panel other than the control point are not generally reliable. Sur-face and near-surface values of velocity away from the control points are usually obtained by interpolation.

Parabolically curved panels have been used in at least three higher-order, subsonic-flow codes (refs. 31-33). For a given number of panels on a curved surface, such panels generally provide better accuracy than do flat panels. The parabolically curved panels are similar to the flat panels in that they still leave gaps between adjacent panels.

As indicated in the introduction, the difficulties encountered with extending panel methods to supersonic flow were eventually traced to using doublet-strength dis-tributions that were not continuous at panel edges. The task of building continuous donblet-strength distributions requires that the panel edges themselves be contiguous (have no gaps) since the doublet strengths are referenced to the panel geometry.

There are several ways to connect an MxN array of arbitrary points with (M-1)x(N-1) quadrilateral panels that have no gaps along adjacent panel edges. Probably the simplest is to use several piecewise flat subpanels to connect four corner points. The most recent version of the QUADPAN code uses four triangular subpanels (fig. 12(a)), whereas the PAN AIR, HISSS, AND MCAERO codes use four triangular subpanels and a flat interior parallelogram (fig. 12(b)). Although not shown in figure 12(b), the flat parallelogram portion of the panel is further subdivided into four more triangular panels defined by the diagonals of the parallelogram. This is done to facilitate the construction of the continuous quadratic doublet-singularity distributions within the panel.

Two codes, SOUSSA and the most recent version of the LEV code (leading-edge vortex code, also called the flee vortex sheet code, ref. 27), use hyperbolic-paraboloid panels, (fig. 12(c)). The curved shape of these panels is the same as for a structural plate that is loaded as in fig-ure 12(d), that is, a twisted shape is produced (ref. 34). The side edges remain straight and hence produce no gaps.

Use of the contiguous piecewise fiat panels has enabled the influence-coefficient-equation surface inte-grals (eqs. (9) and (11)) to be evaluated analytically for both subsonic and supersonic free-stream flow. Use of the hyperbolic-paraboloid panels has been restricted to sub-sonic flow (no one has been able to evaluate the surface integrals when using these panels for supersonic flow).

Singularity Distributions

Most codes use a combination of source and doublet distributions on the panels. The primary exceptions are the Woodward-Carmichael and USSAERO (also called Woodward-2) codes, which use elementary horseshoe vortices instead of doublets. The elementary horseshoe vortex potential is obtained by integrating the doublet potential in the x (streamwise) direction (ref. 35, p. 87). The strength of the vortex distribution, also called the sur-face vorticity vector Y, is related to the doublet strength as follows (ref. 2, eq. (A-3); ref. 11, eq. (03.3.9)) (eq 18). Both \arrow{\gamma} and \arrow{\delta}\mu are in the plane of the panel, \hat{n} is a unit vector normal to the panel, and \cross is the cross-product operator.

Source, doublet, and vorticity surface distributions cause certain flow properties to be discontinuous with respect to the panel surface; that is, there is a jump in the flow properties. In general, this means that the velocity vectors on opposite sides of a panel are different. This is illustrated in figure 13, where the subscripts U and L refer to the "upper" and "lower" sides of the panel, respectively. For incompressible flow, a source-only panel causes the normal component of the velocity at point (\zeta, \eta) of the panel to jump by an amount given by (eq 19) and a doublet-only panel causes both the potential and the tangential component of velocity to jump. These jumps are given by (eq 20) (eq 21). The jump property given by equation (21) is derived in appendix B. Each of the jump properties is derived in reference 33.

The velocity jump across a panel for incompressible flow is then given by (eq 22) The generalization of equation (22) for compressible flow is given by equation B.3.29 of reference 12 (version 3.0, 1990).

Equations (18) and (22) provide some guidance for selecting consistent functional forms for the singularity distributions. A constant-strength source panel combined with a constant-strength vorticity panel or a linear-strength doublet panel will produce a velocity-jump dis-tribution that is constant over the panel surface. A linear velocity-jump distribution requires a linear-source and linear-vorticity or quadratic-doublet distribution. As will be discussed later, popular sets of boundary conditions are those that include setting the perturbation potential to zero on the interior side of panels (the side not wetted by the external flow field). Consequently there is no interior perturbation velocity, and the velocity jump across the panel equals the perturbation velocity on the exterior side.

As discussed in reference 32, consistency also depends on the surface shape used for the panels. Flat panels are consistent with constant sources, and parabolic panels are consistent with linear sources. In this sense, the methods of reference 32 and 33 (which use parabolic panels, linear sources, and linear vorticity or quadratic doublets) are among the few totally consistent higher-order formulations.

Continuity of Doublet Strength

The reason discontinuous doublet strength can cause disastrous numerical problems in supersonic flow is illus-trated in figure 14. For a problem simple enough to solve in closed form, the actual doublet strength might look something like that in the top part of the figure, that is, \mu varies in a continuous manner. A panel method that does not enforce doublet-strength continuity at panel edges will produce a solution like that shown in the lower part of the figure, that is, the doublet strength jumps at the panel edge. Now consider either of the panels separately. In this view the doublet strength has a nonzero value at the panel edge. The velocity field produced by the doublet distribu-tion on either panel is given by the second term of equa-tion (11). This term can be integrated by parts (appen-dix B of ref. 11) to obtain (eq 23) where the line integral is around the panel perimeter.

The surface integral in equation (23) involves the sur-face vorticity vector _, and is called the regular part of the doublet velocity. The line integral involves the doublet strength at the panel edges and is called the line-vortex term. The two line-vortex terms from the common edges of the two panels in figure 14 produce a single line-vortex of strength F(/) = A_t(t), where AIx(l) is the jump in the doublet strength along the panel edge. The velocity field from this line vortex is spurious since the actual doublet strength does not jump at the panel edge; that is, the dis-cretization has introduced velocities that should not exist.

In subsonic flow these spurious velocities decay rapidly with distance from the edge and usually do not cause serious problems. In supersonic flow, these veloci-ties persist, their effect propagating down the Mach cones. Consequently, erroneous incremental flows continue to exist at control points within the domain of influence of the disturbance points, thereby introducing errors in the AIC matrix. These errors are frequently serious enough to produce a totally incorrect solution for the flow.

Now, if it were known in advance that the discretized doublet strength would always be continuous in value, then the resultant line-vortex term along common panel edges of adjacent panels would always be zero. Hence, the line-vortex term in equation (23) would not even have to be computed, but could be simply discarded. This approach is taken in the subsonic-supersonic-flow PANAIR and HISSS codes: the doublet strengths are made continuous by construction and the line-vortex terms are thrown away.

A comparison of results for a 3%-thick swept wing at M** = 2.05, based on discontinuous and continuous quadratic-doublet distributions, is shown in figure 15 (figs. 9 and 10 in tel 36; figs. 48 and 62 in ref. 37). In the discontinuous-distribution case, there is a pressure spike at about 80% chord of the station just inboard of the tip. This is caused by a doublet-strength discontinuity near the leading-edge region of the wing tip, the effect of which propagates down the Mach cone originating at the tip leading edge. The spike is not present in the continuous doublet-distribution case.

Summary of Methods

Table 1 summarizes the basic features of several panel codes. Codes that handle only subsonic flow are listed separately from those that treat both the subsonic and supersonic cases. The QUADPAN code is listed twice because it was originally a subsonic-only code and was later revised to do supersonic flow. The year given next to each code's name is the approximate year the code was introduced. The panel geometry and the singularity types and their spatial distributions are listed for each code. A footnote indicates that geometry or singularity type is con-tinuous from panel to panel. The table includes what are probably the best known, or most generally available codes; other codes (often proprietary) also exist. Refer-ences 39 through 41, as well as other sources (material distributed at the 1985 AIAA Workshop on Aerodynamic Analysis Using Panel Methods), compare results from several of the codes for various configurations.

Codes for subsonic-only flow-- The lower-order Hess code (also known as the Douglas-Neumann code) is considered by many to be the first practical implmentation of the panel method for quite general geometry (ref. 2). It actually exists in several versions, the first being a source-only version that did not treat lifting problems. The lifting case was later added. The Hess family of codes is one of the few that enforce the Kutta condition with a nonlinear pressure rule. As shown in reference 40 this can sometimes give a better trailing-edge pressure result than is possible with a linear approximation to the pressure rule.

MCAERO uses higher-order continuous doublet dis-tributions on piecewise flat continuous panels (tel 17). One of its important capabilities is the use of analytically differentiated (with respect to panel coordinates) influence-coefficients. These are used to efficiently implement the so-called design problem, that is, the prob-lem of determining the geometric shape required to pro-duce a specified pressure distribution.

SOUSSA stands for steady, oscillatory, and unsteady subsonic and supersonic aerodynamics (tel 8). Although it produced some supersonic results by using a small but nonzero reduced frequency, it turned out to be incapable of doing steady supersonic flow. One of its legacies was to popularize the so-called Morino boundary-condition formulation, which is discussed in section 5. This formu-lation led to renewed interest in lower-order subsonic codes and resulted in the development of the VSAERO and QUADPAN codes.

VSAERO (from vortex separation aerodynamics analysis) is one of the few codes that contain a procedure for calculating the shape and location of the trailing-wake system (ref. 18).

The LEV code (leading-edge vortex; ref. 27) is also known as the free-vortex sheet code. It was designed specifically to model the vorticity shed from sharp leading edges of swept wings (fig. 5).

The higher-order Hess code (ref. 32) is the most recent version of the Hess family of codes. It uses parabolically curved panels in conjunction with linear source and linear vorticity (quadratic doublet) distributions.

Codes for subsonic or supersonic flow-- The Woodward-Carmichael code (refs. 4, 5) (also known as the constant-pressure panel code and the Woodward-I code) was described in the introduction. It is still used for simple configurations that can be approximated with the mean-surface representation of figure 1(b).

USSAERO (unified subsonic/supersonic aerodynam-ics) is also known as the Woodward-II code (ref. 7). It differs from the Woodward-I code in that the line distri-butions used for bodies were changed to constant-strength source panels, and the constant-strength vorticity panels were changed to distributions wherein the strength varies linearly in the chordwise direction and is constant in the spanwise direction. The corresponding doublet strength (quadratic-linear) is not continuous, so flat mean-surface models of lifting surfaces are usually required for super-sonic flow.

PANAIR (panel aerodynamics) is generally consid-ered to be the first actual-surface paneling code with reli-able numerics for supersonic flow (refs. 9-14). In addition to using continuous doublet distributions, it was also nec-essary to incorporate forward-weighted splines for the supersonic case. This is somewhat analogous to the use of forward-differencing in the finite-difference computa-tional fluid dynamics (CFD) codes. The basic technology in PAN AIR is also used in the HISSS code (higher-order subsonic/supersonic singularity method, refs. 16 and 42).

QUADPAN (quadrilateral panel aerodynamics pro-gram) started out as a subsonic-only code and used constant-strength sources and doublets (ref. 19). Later, the doublets were changed to a linear-continuous distribution so that supersonic flow could be handled (from 1985 AIAA Workshop on Aerodynamic Analysis Using Panel Methods).

It has been well established that continuous doublet distributions are essential for reliable numerics in super-sonic flow. However, it is probably fair to say that there is not a unanimous agreement between panel code develop-ers about the need for higher-order approaches for the subsonic-flow codes. The advantages that lower-order codes offer over higher-order codes are (1) less work to derive the influence-coefficient equations; (2) simpler coding implementation (because the higher-order approaches must relate information involving different panels, which leads to special cases and logic); and (3) far fewer arithmetic operations, hence lower run costs. One reason some subsonic-only codes use higher-order approaches is that their developers believe the numerics are more reliable than those of the lower-order approaches for highly complex geometry. In the case of two codes that use parabolic panels with linear sources and quadratic doublets, numerical calculations have demonstrated that fewer panels are required than for lower-order methods to obtain a given accuracy (refs, 32, 33). This seems to be especially true for internal flows.

MODELING

Modeling refers to particular techniques used to simulate flow about an object. The modeling tools at the disposal of a panel-method user are (1) the geometric generality that panel codes provide, (2) the use of sources and doublets, individually or in combination, and (3) boundary conditions. In this section we will discuss boundary conditions, their interplay with the source-doublet jump properties introduced in section 4, and how the combination can be used to model in different ways. There is also a brief look at wake modeling.

In section 2 we discussed the general limitations inherent in the Prandtl-Glauert equation derivation, namely, the fl0w is represented as being inviscid, irrota-tional (potential), and linear. If the free-stream Mach number is supersonic, there are additional geometric restrictions; these are illustrated in figure 16.

In supersonic flow," the Prandtl-Glauert equation admits solutions for solid surfaces only if the surfaces are swept back more sharply than the Mach cone. Thus, the higher the Mach number, the more streamlined must be the aircraft. This restriction means that forebodies must be pointed, not blunt. Wing leading edges can be blunt only if they are swept behind the Mach cone (a so-called sub-sonic leading edge).

Surfaces can be swept at angles smaller than the Mach cone only if they do not represent solid surfaces. These so-called superinclined panels (ref. 11) are used for nacelle inlet faces and nozzle exit planes. Superinclined panels use both sources and doublets, and both boundary conditions must be prescribed on the downstream side of the panel. Numerical experiments indicate that these downstream boundary conditions must also be of a spe-cific type, namely, that the potential and its normal deriva-tive must be specified. This is analogous to an initial value problem. Superinclined panels have no upstream influence. The panels generate a downstream flow only and simply absorb any flows that run into them. Refer-ence 11, appendixes A and B in reference 12, and refer-ence 37 are recommended sources of information on the supersonic aspects of modeling.

Mass-Flux and Velocity Boundary Conditions

The physical description of a real flow at a surface is given by the no-slip boundary condition (eq 24a) or (eq 24b) These equations state that the total velocity vector \arrow{V}, or the mass flux vector (\rho\arrow{V}) is zero at a solid surface. For inviscid flow, the tangential component of the velocity cannot be prescribed (unless the pressure is known) and equations (24) are replaced with the zero normal-flow boundary conditions, (eq 25a) or (eq 25b) which must be supplemented with a Kutta condition for subsonic trailing edges.

Equation (25b) is nonlinear since the density is a function of the unknown velocity; a linear approximation to equation (25b) is (from sec 1.11 of ref. 30) (eq 26) Equation (26) involves neglecting terms that are of the same order of magnitude as those neglected in deriving the Prandtl-Glauert equation. Dividing by \rho_{\infty} gives (eq 27a) where, from equation (8) (eq 28) which is linear in the components of the pertubation velocity, and depends on the free-stream Mach number. If the panel is not meant to be impermeable, for example, at the inlet face, then the right-hand side of the equation (28) is replaced with a specified nonzero value.

For a fixed subsonic Mach number, a solution for \phi(x,y,z) can be obtained in either of two ways. The first is to solve the "true" compressible problem, that is, solve the Prandtl-Glauert equation as it stands, using the mass-flux boundary condition applied to the true geometry. The second (and better known) way is to solve the equivalent incompressible problem by using the Prandtl-Glauert transformation to convert equation (1) to Laplace's equation (sec. 7-1 in ref 35): (eq 29) using the velocity boundary condition (eq 30) applied to the transformed geometry. In equations (29) and (30) the primes indicate that all the variables are in the transformed coordinates, for example, ____. Transforming the solution for \phi' back to the true (physical) variables gives the solution for \phi for the true geometry. The important point is that the solution for \phi obtained from the first approach will be the same as that obtained from the second approach (ref. 43). Thus, solving for the flow about the true geometry by using the subsonic Prandtl-Glauert equation with mass flux (instead of velocity) boundary conditions, is mathe-matically the same as solving the equivalent incompress-ible problem with velocity boundary conditions applied to the transformed geometry. The first approach has the advantage that it can also be used for supersonic flow (where there is no equivalent incompressible problem).

Note that the linearized mass flux defined by the equations (27) is not parallel to the true mass flux \rho\arrow{V}, and hence, is not parallel to the velocity vector \arrow{V}=\arrow{V}_{\infty}+\arrow{v}. This inconsistency is one of the prices paid for the linearization. The consequences of this will be illustrated with the following two examples. (Theoretical discussions of mass-flux and velocity boundary conditions appear in refs 44-46).

The first example demonstrates the accuracy of the mass-flow boundary condition and the Prandtl-Glauert equation for supersonic flow over a thin wedge. It also demonstrates the jump properties of the mass-flux and velocity vectors across shocks predicted with the Prandd-Glauert equation.

Panel-method results obtained from the PAN AIR code are shown in figure 17 (from ref. 9). The boundary condition applied to the exterior side of the panels is the linearized max flux condition \arrow{W}\dot\hat{n}=0. This forces the linearized mass flux vector \arrow{W}, instead of the resultant velocity vector \arrow{V}, to be parallel to the wedge faces. Pressures have been computed with the isentropic and second-order pressure rules (eqs. (8.10) and (8.11), respectively, in ref. 23). The PAN AIR results are very close to those of nonlinear shock-expansion theory, having a greater pressure magnitude on the frontward-facing compression surface than on the rearward-facing expansion surface. This contrasts with the classic linear thin-airfoil theory, which predicts equal and opposite pressures on the two inclined faces.

How does one explain that PAN AIR, which solves the linear Prandtl-Glauert equation, produces results that agree more closely with nonlinear shock-expansion theory than with classic linear thin-airfoil theory? Apparently, for this particular problem, the velocities are small enough to make negligible the nonlinear terms in the differential equation. This conclusion is justified by the close agree-ment between the isentropic and second-order pressure rules. (It is good practice to always compute the pressures by these two rules; a substantial disagreement in the results is a sure sign that the small-perturbation assump-tions of the Prandtl-Glauert equation are being violated.) Since classic thin-airfoil theory is also based on the Prandtl-Glauert equation, but in addition depends on complete linearization of the boundary conditions, the explanation must reside in the boundary conditions. Although PAN AIR uses the linearized mass-flow W in the boundary condition fie. fi = 0, the unit normals fi act at the actual wedge surfaces. In the classic linear-theory solution the boundary conditions are not applied to the true geometry. Instead, the boundary conditions are applied along the straight upper side of the wedge, and the flow is made to turn through the wedge angles, as in the mean-surface modeling of the Woodward-Carmichael code (see sec. 1). So, this appears to be an example of an instance in which the solution is more affected by approx-imations made to the boundary conditions than it is by those made to the differential equation.

The mass-flux and velocity vectors at two points on opposite sides of the Mach line emanating from the wedge apex are shown in figure 18. We first note that for expan-sive flow about a corner, such as at the wedge apex, the Prandti-Mayer solution predicts an expansion fan. That is, away from the apex the flow properties change gradually through the fan. However, the Prandtl-Glauert equation predicts abrupt changes that occur along the apex Mach line. This "expansion shock" approximation to an expan-sion fan is a poor representation at large distances from the apex, but is accurate near the apex. Our point here is to illustrate what happens to a panel-method solution across such a shock (whether expansive or compressive). Both the mass-flux and velocity vectors are discontinuous across the Mach line; they jump in both magnitude and direction. For the mass flux, the tangential component jumps, and the normal component is continuous. For the velocity, the tangential component is continuous, and the normal component jumps. This is the same behavior as predicted by the nonlinear Rankine-Hugoniot relations, thus providing some measure of confidence in the lin-earized mass-flux approach.

The second example (fig. 19) shows the pressure at any point on 10\deg and 15\deg half-angle cones at zero angle of attack, for free-stream Mach numbers between 1.0 and 4.0 (fig. 4.1 in ref. 15). At the lower Mach numbers, the mass-flux-boundary-condition solution agrees more closely with the exact Euler-equation solution than does the velocity-boundary-condition solution. As the Mach num-ber increases, the mass-flux and velocity-boundary-condition solutions become less accurate. The mass-flux-boundary-condition solution rapidly diverges from the exact solution, and crosses the less rapidly diverging velocity-boundary-condition solution. As the cone angle increases, the Mach number at which the mass-flux and velocity-boundary-condition solutions cross one another becomes smaller. This supersonic "thick-body" behavior is responsible for the fact that when panel methods are applied to nonslender fighter forebodies, the velocity-boundary condition often gives superior answers to the mass-flux boundary condition in the region just ahead of the canopy (refs. 47, 48). For such cases, the mass-flux boundary condition can actually produce negative pres-sure coefficients, as suggested by the 15 ° cone solution behavior of figure 19. For bodies that are adequately slen-der, the linear Prandtl-Glauert equation with mass-flux boundary conditions can provide good answers, for exam-ple, the Mach 1.6 pressures on the B-1 forebody presented in references 36 and 37. The important conclusion to be drawn from figure 19, and from other examples in ref. 15, is that inviscid supersonic flow solutions based on the Prandtl-Glauert equation can be substantially in error if the Mach number, thickness, or angle of attack are too large. In such cases, codes based on nonlinear theories (for example, the full-potential equation or Euler equa-tions) must be used for reliable answers.

Interior Potential

An actual-surface panel model of an aircraft generally produces a set of panels that separates space into two or more distinct regions: enclosed interior volumes and an external volume extending from infinity to the external side of the panels. The flow in the external volume corre-sponds to the physical flow field being modeled. The flow in the interior volumes is fictitious but, as will be shown, can be used to advantage.

Newcomers to panel methods often find the idea of an internal flow field to be strange in that no such flow exists inside a real wing. However, it must be remembered that we are using sources, doublets, and boundary conditions to create the flow fields, and flow will exist on either side of the source-doublet panels. The external-internal flow fields are in general independent of one another; they depend on the boundary conditions on the external-internal sides of the panels. These include direct boundary conditions explicitly imposed on each side of the panel and indirect boundary conditions resulting from some specification of the source or doublet strength. The indirect boundary conditions arise because of the jump in flow properties produced by sources and doublets.

The boundary condition implied by the equation (20) jump property is the cause of a common mistake made by new users of panel methods, namely, not specifying the potential at a point of an interior volume. When an interior volume is totally enclosed by panels that all have doublet distributions, the potential must be specified at one or more interior points, and this must be done with a source distribution. If this is not done, the problem is ill-posed, that is, it does not have a unique solution, and the AIC matrix will theoretically be singular (will not possess an inverse). (In practice, numerical round-off error often changes a theoretically singular matrix to just an ill-conditioned matrix.)

The reason for the problem being ill-posed is generally explained as follows. At \arrow{x}=\infty the value of the potential owing to a source or doublet is zero (this is called the far-field boundary condition). At the exterior side of panels enclosing an interior volume the potential will have some distribution \phi_{ext}(x,y,z) that depends on the surface boundary conditions and on the resulting source-doublet strength distribution. On the interior side of the panels, the potential is \phi_{int}(x,y,z), which by equation (20), differs from the external distribution by an amount equal to the doublet-strength distribution \mu(\zeta,\eta). Now, if at some interior point, a constant c, arbitrary except that it must satisfy any boundary condition associated with the point, were to be added to \phi_{int}(x,y,z), then \phi_{int}(x,y,z)+c would also satisfy the Prandtl-Glauret equation and, hence, be another solution. The way to make the interior solution unique is to specify a value for the potential at an interior point, thus determining the constant c. (The same argument can be made for the exterior side of the panels, which requires that the arbitrary constant be zero to satisfy the far-field boundary condition \phi_{\infty}=0. Thus, no explicit specification of the potential at an exterior point has to be made.)

The interior potential at a point must be specified with a source panel, and not a doublet panel, because of the "physically" different behavior between sources (or sinks) and doublets. With any paneled geometry, the boundary conditions are enforced only at the control points. Hence, when a finite number of control points is used to impose boundary conditions there will be some degree of "leakage." That is, the integral of the source strength (jump in normal component of mass flux) over the surface enclosing the interior volume will not be exactly zero. When a source panel is used to specify velocity potential in the interior domain, the panel is capable of generating or consuming fluid to conserve mass that is leaked out of or into the interior domain. (The source strength is one of the unknowns that are solved for.) The doublet panel does not have the ability to generate or swallow fluid (at least, in a net sense), and consequently is not capable of handling the leaked flow. Numerical examples demonstrating this behavior are given in cases 5 and 6 of reference 49 (pp. 6-7).

Another interesting example in reference 49 is case 7. Here, a thick wing is modeled using doublets-only for the wing upper and lower surfaces and sources-only for the panels closing off the wing-tip opening. Using only a single type of distribution on all panels allows only one boundary condition per panel, which was chosen to be \arrow{W}\dot\hat{n}=0 on the exterior panel sides. The interior potential was not specified anywhere inside the wing. At first glance this also apperas to be an ill-posed problem, but it is not. The source-only tip-closure panels do not produce a potential jump. Therefore, as far as the potential is concerned, this model behaves as if there were no tip closure, and in this sense the wing interior is not a separate domain from the exterior domain. Hence, the constant \phi_{\infty}=0 applies to both the exterior and interior part of the single domain. The zero-normal-flow boundary condition applied at the tip, however, does separate the flow into two regions.

In actual practice, the interior potential is set not at just a single point, but at the interior-side control point locations of all panels. If the boundary condition \phi_{int}(x,y,z)=0 is used at all such control points there will be no internal pertubation flow (assuming an adequate number of panels), giving uniform free-stream flow in the interior volume. If both sources and doublets are used with each panel, then another boundary condition is needed, and it is selected to control the flow on the exterior sides of the panels. Two alternative ways of doing this, for flow about a solid (impermeable) aircraft surface are described next. The first is a direct approach, the second is an indirect one.

The direct approach is illustrated in figure 20(a). On the exterior panel side, the zero normal mass flux boundary condition (eq. (28)) is imposed. The resulting boundary condition pair at each panel control point is (eq 31a) (eq 31b) The indirect approach (sometimes referred to as the Morino formulation) is illustrated in figure 20(b). Here, equation (31a) is replaced by a condition imposed o the source strength, and the boundary condition pair is (eq 32a) (eq 32b) The equivalence between these two boundary condition pairs can be shown by starting with the definition of the source strength (see eq. (8)): (eq 33)

If the pertubation potential is zero everywhere in the interior region, then its derivative is zero in every direction and \arrow{w}_{int} is then zero. (Equation (32b), applied to the interior side of every panel normally produces a good approximation to this condition.) Then equation (33) becomes (using the definition for \arrow{W} given by eq. (27a)) (eq 34) Now imposing the \arrow{W}_{U}\dot\hat{n}=0 condition sought (eq. (31a)), we get equation (32a).

The indirect approach has two features that reduce run cost. First, only the influence-coefficients for the potential, and not the three components of velocity, have to be computed. Second, since the source strengths are specified, only the doublet strengths have to be solved for, cutting the size of the [AIC] matrix roughly in half. (The reason the size is not necessarily cut exactly in half is that some codes use more doublet unknowns than source unknowns, e.g., for wakes, and for doublet matching in the continuous-doublet codes.)

The direct approach requires that both potential and velocity influence-coefficients be computed, and both source and doublet strength must be solved for. For a given number of panels this causes the direct approach to be more costly than the indirect approach. The advantage of the direct approach is that it is sometimes more accurate than the indirect approach. This is apparently because the indirect approach depends on the internal pertubation potential being zero everywhere for equations (32) to accurately represent equations (31). In practice, the potential is set to zero only at discrete interior points; consequently, the indirect formulation sometimes produces more error than does the direct approach.

An additional advantage, common to both formula-tions, is the property that the internal flow is everywhere uniform (it equals the free-stream velocity). For super-sonic free-stream Mach numbers this prevents the forma-tion of internal disturbances that otherwise would propa-gate along Mach lines and reflect from the internal sides of the panels. If not eliminated, these reflections can cause severe internal flow disturbances that are "felt through the panels" and degrade the external flow-field solution. Although there are no Mach-line disturbances in subsonic flow, large velocities are normally produced by the line-vortex behavior of constant-strength doublet panels. These large velocities are eliminated in the internal flow region by the internal potential being set to zero. This is apparently why solutions obtained with the newer constant-strength source-doublet formulations are not as sensitive to the panel layout as were earlier constant-strength codes (which did not use a zero interior-perturbation potential).

A final example uses doublet-only panels to represent a wing, including the tip closure. Since the interior volume is totally enclosed with doublet panels, the potential must be set in the interior volume for the problem to be well-posed. Also, having only one type of singularity dis-tribution means that only one boundary condition can be employed per control point. The direct approach would be to set the interior potential with an additional source panel located inside the wing, and specify \arrow{W}\dot\hat{n}=0 on the exte-rior. A clever indirect approach is to simply set the total (instead of the perturbation) potential to zero on the inte-rior side of each panel. For sufficiently dense paneling, this will make the total potential close to uniformly zero at all interior points. Hence the gradient of the total potential in any direction, which is the total (linearized) mass-flux in that direction, will also be zero. Consequently, the total mass-flux component normal to the interior side of the panels will be zero. Then, since doublets do not produce a jump in normal mass flux, the normal mass flux on the exterior sides of the panels will also be zero, that is, \arrow{W}\dot\hat{n}=0 on the exterior is produced.

Wakes

Wake panels are used to enforce the Kutta condition at sharp (usually the trailing) edges of lifting surfaces. Since wakes trail downstream from these edges, they also influence the flow experienced by downstream components. For example, the load experienced by a wing will depend on the proximity of the wake from an upstream canard.

Figure 21 shows a PAN AIR model of a supersonic fighter wing-canard concept (ref. 50). Two of three wake models tried for this configuration are shown in figure 22. Simple, flat wakes for both the wing and canard are used in the first model, the wakes being positioned approximately in the wing-canard mean plane (fig. 22(a)). For each angle of attack, the resulting span loadings were then used by a multiple lifting-surface code to compute the approximate rolled-up shape of the canard wake. The rolled-up canard wakes were used for the second model (fig. 22(b)). The third model (not shown) again used flat wakes, but for each angle of attack the canard wake was aligned with the free-stream velocity direction. The results in reference 50 indicate that for this configuration, the second and third models give essentially the same lift and moment, and these agree more closely with wind-tunnel results than do the results from the first wake model. Lift and pitching moment coefficients at Mach 2.2, obtained from the third wake model, are shown in figure 23. There is a noticeable difference in the PAN AIR-predicted lift-curve slopes for the isentropic and second-order pressure rules, indicating that nonlinear effects are beginning to become important. Also shown are results from the USSAERO (version B) program; it does less well at predicting the moment.

Another wing-canard configuration is shown in fig-ure 24 (ref. 51). The effect of aligning a fiat-canard wake model with the free stream instead of with the chord plane (unaligned position of fig. 25) for this close-coupled wing-canard is illustrated in figure 26. This figure shows the spanwise circulation distribution of the canard and the wing for the aligned and unaligned fiat-canard wakes. The canard-wake position does not significantly affect the lift distribution along the canard. It does, however, have an important effect on the wing-lift distribution. Moving the canard wake from the unaligned position to the aligned position causes two major changes in the flow over the wing, the dominant one of which is an increase in wing lift inboard of the canard-tip station. This is due to the diminished canard-wake downwash field, raising the effective angle of attack of the inboard wing section. The secondary effect is a loss in wing lift outboard of the canard-tip station. This loss in lift is due to diminished spanwise velocity imparted on the upper surface of the wing by the canard wake, which is due to the increased distance between the canard wake and the wing. The net effect of these changes caused by the aligned wake model is an increase in total lift over that predicted by the unaligned wake model.

The final example, taken from reference 52, is a PAN AIR flaps-down analysis of the Boeing 737-300 (fig. 27). The actual and computational flap geometries for the flaps-15 setting, used for most takeoffs, are shown in figure 28(a). The assumed flap, wing, and slat wake positions are shown in figure 28(b). Leading- and trailing-edge spanwise geometry discontinuities owing to the deflected flaps and slats also had to be treated; these are discussed in reference 52.

The drag buildup was obtained from a combination of methods. The profile drag for everything but the wing was estimated from handbook data. The wing profile drag was estimated from a two-dimensional multielement-airfoil panel code coupled with a viscous model (ref. 53). The induced drag was obtained from PAN AIR's surface-pressure integrations. Figure 29(a) shows the drag contributions for both the flaps-15 (F15) setting, and a lower setting called flaps-1 (F1). The resulting computational lift-to-drag ratios are compared with flight data in figure 29(b).

TRANSONIC FLOW

Panel-method codes have existed for about 25 years and are still the only codes routinely used to analyze flow about complex three-dimensional configurations. Their general inability to solve nonlinear problems is a serious drawback, however, and most current CFD research involves finite-difference, finite-volume, or finite-element approaches to solving nonlinear flow equations. In the United States, the finite-difference and finite-volume approaches currently seem to be in most favor. Unfortu-nately, most finite-difference and finite-volume approaches require well-structured flow-field grids that conform to the surface of an aircraft. It "is the difficulty in generating 'suitable grids'" (ref. 54) that is the major technical obstacle to routinely computing inviscid tran-sonic flow about realistic aircraft. In contrast, it is rela-tively easy to produce panel-method-type grids that are only on the aircraft surface.

This raises the question: Is it possible to combine the surface geometry grid used by panel methods, with some easy-to-generate flow-field grid and use the combination to solve nonlinear fluid flow problems? The answer is yes, as demonstrated by the TRANAIR code and the work that led to it (refs. 55-59). This section provides a brief introduction to the technique, and how it evolved. To the user, TRANAIR appears to be a panel code since the input is panel-code-like; however, the solution techniques TRANAIR uses are not those of a panel code.

The basic approach used in TRANAIR is to embed the surface panels in a rectangular box of grid points, as shown in figure 30. The initial formulation (sec. II of ref. 56) was based on Green's third theorem (eq. (7)), and combined the surface-integral-generated influence-coeffi-cients with the volume integrals. The rectangular grid was used to evaluate the volume integral, for every point of the rectangular and surface grids, with fast Fourier trans-forms. This Green's theorem approach was able to solve the nonlinear fuU-potential equation if there were no shocks, but proved to be unstable when shocks were pres-ent. By resorting to more powerful mathematical methods involving optimization, supercritical results using influence-coefficients were subsequently obtained (refs. 55, 56); however, the computational cost was extremely high.

The most recent version of TRANAIR (refs. 57-59) still uses panels and the rectangular grid, but does not use influence-coefficients. Instead, the cells formed by the rectangular grid are used to discretize the full-potential equation with tri-linear basis-function finite elements. The surface panels, which slice through some of the finite-element cells, alter the finite element discretization in teh vicinity of boundary surfaces.

At the perimeter of the rectangular grid, the equation set being solved changes from the full-potential equation to the Prandtl-Glauert equation. As a consequence, the rectangular computational grid need only encompass the nonlinear flow regions (which are only near the aircraft). The far-field boundary condition of zero potential at infin-ity is automatically satisfied by the discrete Green's func-tion (for the Prandtl-Glauert equation) used in the formu-lation. Consequently, the solution domain extends to infinity, even though the computational grid is finite, as indicated in figure 30.

The finite-element discretization yields a set of non-linear algebraic equations which are solved iteratively, using Newton linearization, multiple preconditioners, and an optimzation algorithm called GMRES (generalized minimal residual). Details are given in references 58-60. The computer run cost based on this approach is less than with the original influence-coefficient/optimization approach.

The input to TRANAIR is essentially that of the PAN AIR code, that is, the surface grids of panel corner points supplemented by the box of rectangular flow-field grid points. This enables transonic flow to be computed about very complex configurations without having to generate a surface-conforming flow-field grid. For example, TRANAIR has been used to compute transoinc full-potential solutions for the F-16A, using the rectangular grid box and paneled geometry of figure 30. For the half geometry shown, the grid box contains 129x33x33 points, and the aircraft contains about 3,500 panels. Two views of the surface paneling are shown in figure 31. In this model, the wing-tip missles and launchers are not included.

F-16A supercritical wing-pressure results (ref. 57) are shown in figure 32. The free-stream Mach number is 0.9 and the angle of attack is 4\deg. The experimental data for the outboard station indicates separated flow near the trailing edge. This is probably a result of the wing-tip missiles and launcher that were part of the wind-tunnel model. The wind-tunnel data indicate a shock at approximately 75% chord for the four inboard stations. TRANAIR also indi-cates a shock, but it is slightly downstream of the shock predicted by the wind tunnel. This result is generally expected from a conservative full-potential solution, a result of the absence of a boundary-layer correction (fig. 40 in ref. 61). The shock predicted by TRANAIR is smeared over 5 to 6 grid cells, and at the 2 outboard sta-tions where there are only about 11 rectangular-grid points spanning the chord, the shock is completely washed out. Increasing the grid density in the x-direction (streamwise) greatly reduces the shock smearing. To bet-ter resolve rapidly varying flow behavior without increas-ing the flow-field grid density everywhere, techniques such as local grid refinement, wherein individual cells of the rectangular grid are subdivided into smaller cells, or higher-order finite-element-basis functions are needed.

The coarseness of the rectangular grid is particularly evident in the leading-edge region of the wing. An example is shown in figure 33, at about 70% semispan. The first 18% of the wing chord is spanned by only three cells of the rectangular grid. Consequently, with this coarse uniform grid, TRANAIR fails to capture the leading-edge pressure peaks. The subsequent addition of local grid refinement has enabled these pressure peaks to be resolved (ref. 62).

Panel Methods -- An Introduction

Summary