Using An Adjoint Approach to Eliminate Mesh Sensitivities in Computational Design Eric J. Nielsen* and Michael A. Park† NASA Langley Research Center, Hampton, Virginia, 23681 An algorithm for efficiently incorporating the effects of mesh sensitivities in a computa- tional design framework is introduced. The method is based on an adjoint approach and eliminates the need for explicit linearizations of the mesh movement scheme with respect to the geometric parameterization variables, an expense that has hindered practical large-scale design optimization using discrete adjoint methods. The effects of the mesh sensitivities can be accounted for through the solution of an adjoint problem equivalent in cost to a single mesh movement computation, followed by an explicit matrix–vector product scaling with the number of design variables and the resolution of the parameterized surface grid. The accu- racy of the implementation is established and dramatic computational savings obtained using the new approach are demonstrated using several test cases. Sample design optimizations are also shown. Nomenclature D = vector of design variables f = cost function K = linear elasticity coefficient matrix L = Lagrangian function L⁄D = Lift-to-drag ratio Q = vector of dependent variables R = discretized residual vector X = computational mesh L = vector of adjoint variables I. Introduction I N recent years a concerted effort has been made to bring higher fidelity, physics-based computational fluid dy- namics (CFD) simulations into the aircraft design process. Such tools have typically been used to perform validations of designs derived through the use of lower fidelity methods, or in some cases, used in heuristic design methods that rely heavily on experience. However, these advanced CFD tools are now routinely targeted as primary components of automated optimization frameworks. In the field of gradient-based design, one challenge in utilizing solvers based on the Euler or Reynolds-averaged Navier-Stokes equations has been an efficient and accurate method for computing the sensitivity derivatives required by many optimization algorithms. Approaches such as finite differencing,1-3 direct differentiation,4-10 and the com- plex variable method11-16 can be used for calculating these derivatives; however, their cost scales directly with the number of design variables. For typical aerodynamic design problems where this value may be on the order of tens to hundreds, this limitation precludes the use of such methods. To alleviate the computational burden associated with problems containing many design variables, recent work has focused on the use of adjoint methods.17-35 Adjoint methods may be implemented in either a continuous or discrete context, depending on the order in which the differentiation and discretization operations are performed. The relative merits of the two approaches are the subject of much debate in the literature; one such difference is the focus of the current work. Both the continuous and discrete adjoint approaches introduce an auxiliary, or adjoint, variable that is determined through the solution of an additional linear system of equations. Using the same solution procedure as is used for the governing equations, asymptotic convergence rates of both systems are similar, as the eigenvalues of the two systems are the same. Moreover, for the discrete adjoint variant, if the solution algorithm itself is constructed in a manner which is discretely adjoint to the baseline scheme, the asymptotic convergence rates are guaranteed to be identi- cal.19,20,31 *Research Scientist, Computational Modeling and Simulation Branch, MS 128, Senior Member AIAA. †Research Scientist, Computational Modeling and Simulation Branch, MS 128, Member AIAA. 1 American Institute of Aeronautics and Astronautics The primary advantage of the adjoint approach is that the expense does not scale with the number of design vari- ables. Indeed, the solution of the adjoint system itself is independent of the number of design variables. In practice, however, the discrete approach dictates that the effects of the mesh be accounted for in the sensitivity derivatives, whereas these terms do not appear in the derivation of the continuous operators. While perhaps counterintuitive, this difference is actually viewed as one of the primary advantages of the discrete approach where, by definition, the sen- sitivity analysis is discretely consistent with the flow simulation. If these grid-related contributions are not included, it has been shown that the resulting sensitivities can be extremely inaccurate and even of the incorrect sign.29,30 The continuous approach should converge to the discrete result in the case of a sufficiently refined grid, but this condition is seldom, if ever, met in current practice for realistic three-dimensional complex geometry computations. Unfortunately, the expense of computing these grid-related terms has hindered large-scale application of the dis- crete approach. For a single design variable, the use of a direct mode of differentiation to obtain the derivative of a mesh movement scheme based on linear elasticity was shown in Ref. 30 to cost as much as 30% of a flow simulation. Clearly, the expense associated with performing this operation for several dozen design variables will dominate the overall computation. In the current work, an adjoint formulation is introduced that obviates the need for an explicit computation of the grid sensitivities. Such an approach has previously been used in the automatic differentiation (AD) community;36 however, no explicit formulations for handcoding such a technique have been found in the literature. The scheme pre- sented in the current work significantly reduces the expense associated with the mesh linearizations. Demonstrations of the method show that a rigorous discrete sensitivity analysis for problems based on the Navier-Stokes equations may ultimately be performed at a cost comparable to that of the analysis problem. II. Mesh Sensitivities via Forward Mode Differentiation The discrete adjoint technique for sensitivity analysis can be derived in several ways. Here, the approach taken in Ref. 28 is used. Consider the vector of discretized residual equations R for the Euler or Navier-Stokes equations as a function of the design variables D, computational mesh X, and flowfield variables Q. Given a steady-state solution of the form R(D,Q,X) = 0, a Lagrangian function L can be defined as L(D,Q,X,L ) = f(D,Q,X)+L TR(D,Q,X) (1) where f(D,Q,X) represents a cost function to be minimized and L represents a vector of Lagrange multipliers, or adjoint variables. Differentiating this expression with respect to D yields the following: dL (cid:236) ¶ f ¶ X T¶ f(cid:252) ¶ Q T(cid:236) ¶ f ¶ R T (cid:252) (cid:236) ¶ R T ¶ X T ¶ R T(cid:252) = (cid:237) + (cid:253) + (cid:237) + L (cid:253) +(cid:237) + (cid:253) L (2) dD (cid:238) ¶ D ¶ D ¶ X(cid:254) ¶ D (cid:238) ¶ Q ¶ Q (cid:254) (cid:238) ¶ D ¶ D ¶ X (cid:254) Since the vector of adjoint variables is essentially arbitrary, the coefficient multiplying [¶ Q⁄¶ D]T can be eliminated using the following equation: ¶ R T ¶ f L = – (3) ¶ Q ¶ Q Eq. 3 represents the discrete adjoint equation for the flow simulation. The solution of this linear system of equations for three-dimensional turbulent flows on unstructured grids has been demonstrated previously in Refs. 28-31. Once the solution for L has been formed, the remaining terms in Eq. 2 can be evaluated to give the desired sensitivity vec- tor: dL (cid:236) ¶ f ¶ X T¶ f(cid:252) (cid:236) ¶ R T ¶ X T ¶ R T(cid:252) = (cid:237) + (cid:253) +(cid:237) + (cid:253) L (4) dD (cid:238) ¶ D ¶ D ¶ X(cid:254) (cid:238) ¶ D ¶ D ¶ X (cid:254) The [¶ X⁄¶ D] terms in Eq. 4 represent the mesh sensitivities. In Ref. 30, a mesh movement strategy based on the equations of linear elasticity is described. Although generally not as expensive as a flowfield computation, the cost associated with solving these equations is not trivial. If the system is posed as KX = X , (5) surface then the mesh sensitivities may be computed from the following: 2 American Institute of Aeronautics and Astronautics ¶ X (cid:230) ¶ X(cid:246) K¶ D = Ł ¶ Dł (6) surface Note that the solution of this linear system is equivalent in cost to that of the mesh movement (Eq. 5), and must be ob- tained once for each design variable in D. III. Adjoint Approach for Eliminating Mesh Sensitivities Now reconsider Eq. 1, where a subscriptf has been appended to L to indicate the adjoint variable associated with the flow equations, and an additional adjoint variable L multiplying the residual of the grid movement problem has g been introduced: L(D,Q,X,L ,L ) = f(D,Q,X)+L T R(D,Q,X)+L T (KX–X ) (7) f g f g surface Linearizing with respect to D as before yields dL = ¶ f + ¶ R TL + ¶ Q T(cid:237)(cid:236) ¶ f + ¶ R TL (cid:253)(cid:252) + ¶ X T(cid:237)(cid:236) ¶ f + ¶ R TL +L T K(cid:253)(cid:252) –L T ¶ X (8) dD ¶ D ¶ D f ¶ D (cid:238) ¶ Q ¶ Q (cid:254)f ¶ D (cid:238) ¶ X ¶ X f g (cid:254) g ¶ D surface As before, the coefficient multiplying [¶ Q⁄¶ D]T can be eliminated by satisfying the adjoint equation for the flow simulation, Eq. 3. In a similar fashion, the term multiplying [¶ X⁄¶ D]T can also be eliminated by satisfying a second adjoint problem: KTL = –(cid:237)(cid:236) ¶ f + ¶ R TL (cid:253)(cid:252) (9) g (cid:238) ¶ X ¶ X (cid:254)f With the solution of Eqs. 3 and 9, the final form of the sensitivity vector becomes dL = ¶ f +L T ¶ R –L T ¶ X (10) dD ¶ D f ¶ D g ¶ D surface With the formulation outlined above, a single solution of the linear system given by Eq. 9 is required for each func- tion f. If f is some quantity composed of aerodynamic coefficients such as lift, drag, or moments, several observa- tions can also be made about Eq. 10. For design problems in which the shape is held constant and only global param- eters such as the angle of attack or Mach number are allowed to vary, Eq. 9 need not be evaluated, and the third term in Eq. 10 is identically zero. Conversely, for problems involving solely geometric parameterization variables, the first and second terms in Eq. 10 are identically zero, as there is no explicit dependence of f or R on D. In addition, the term L T [¶ X⁄¶ D] is very cheap to compute and only requires an explicit inner product dimensioned by the g surface size of the surface mesh for each design variable. Implementation The software developed in Refs. 28–30 has been extended to include the formulation outlined above. In the previ- ous implementation, the matrix–vector product [¶ R⁄¶ X]TL has been stored to enable rapid computation of each f inner product with ¶ X⁄¶ D, as these residual linearizations are constant for all shape parameters in D. These me- chanics are now used to construct the right-hand side of Eq. 9 in the current work. The parameterization schemes employed here are described in Refs. 37 and 38, and rely on a free-form deforma- tion technique to provide a compact set of design variables for a wide range of configurations. Given the current vec- tor of design variables D, the methods are used to determine the current location of the surface grid points as well as their analytic derivatives with respect to D, [¶ X⁄¶ D] . The schemes are very inexpensive and can be used to surface consistently parameterize families of computational meshes suitable for multiple disciplines. The matrix K is formed using the method of Ref. 30, and is merely transposed once all of the contributions have been included. The system given by Eq. 9 is then solved using the Generalized Minimal Residual (GMRES) algorithm39 implemented in Ref. 30. Since the eigenvalues of K remain unchanged by the transpose operation, the solution for L converges similarly to that of the mesh movement scheme. In the current implementation, adjoint so- g lutions for multiple functions f may be computed simultaneously as outlined in Ref. 31 by storing multiple right- 3 American Institute of Aeronautics and Astronautics 100 10-2 Direct Adjoint 10-4 10-6 al10-8 u d si Re10-10 10-12 10-14 10-16 10-18 Figure 1 ONERA M6 grid used for evaluating 0 20 40 60 80 100 KrylovVector linearizations. Figure 2 Convergence rates for direct and adjoint modes. hand sides for Eq. 9. Once the solution for L has been computed, an explicit inner product with the surface mesh g sensitivities for each design variable yields the final sensitivity vector. IV. Consistency of Linearization To verify the accuracy of the implementation, derivatives of the lift and drag coefficient are computed for fully tur- bulent flow over the ONERA M6 wing geometry40 shown in Fig. 1. The grid for this case contains 16,391 nodes and 90,892 tetrahedra. The freestream Mach number is 0.84, the angle of attack is 3.06 degrees, and the Reynolds number is 1 million based on the mean aerodynamic chord. The computations have been performed in a domain-decomposed environment on 12 Pentium IV processors using three different approaches as outlined in Table 1. All equation sets have been converged to machine accuracy. The first linearization method is a direct mode of differentiation using complex variables. This approach was origi- nally suggested in Refs. 12 and 13, and was first applied to a Navier-Stokes solver in Ref. 11. The primary advantage of this approach is that true second-order accuracy may be obtained by selecting step sizes without concern for sub- tractive cancellation error typically present in real-valued divided differences. Through the use of an automated scripting procedure as outlined in Ref. 41, this capability can be immediately recovered at any time for the baseline code. For computations using this method, the complex step size has been chosen to be 1· 10–50. The second tech- nique used to verify the linearizations relies on the handcoded implementation29,31 of the discrete adjoint system given by Eq. 3 for the flow equations and the handcoded direct differentiation of the mesh terms30 as formulated in Table 1 Schemes used to obtain sensitivities. Method Flowfield Linearization Mesh Linearization Direct differentiation with Direct differentiation with 1 automated complex variables automated complex variables Handcoded discrete adjoint Handcoded direct differentiation 2 (Eq. 3) (Eqs. 4 and 6) Handcoded discrete adjoint Handcoded discrete adjoint 3 (Eq. 3) (Eqs. 9 and 10) 4 American Institute of Aeronautics and Astronautics Table 2 Comparison of sensitivity derivatives for lift and drag coefficients using various approaches. Design Variable Objective Method Function Thickness Shear Camber Twist 1 -0.584383430968430 -0.073891855284066 1.843734584180741 -0.022010251214990 C 2 -0.584383430967291 -0.073891855283698 1.843734584180810 -0.022010251215005 L 3 -0.584383430968115 -0.073891855283921 1.843734584180955 -0.022010251214989 1 0.058894900355748 -0.006835640271421 0.064393773359690 -0.001817294278046 C 2 0.058894900355806 -0.006835640271405 0.064393773359692 -0.001817294278046 D 3 0.058894900355780 -0.006835640271392 0.064393773359720 -0.001817294278046 Eqs. 4 and 6. The third approach is similar to that of Method 2; however, the present adjoint formulation given by Eqs. 9 and 10 is used to evaluate the mesh terms. Sensitivity derivatives for the lift and drag coefficients using several shape parameters located at the midspan of the wing are shown in Table 2. All three methods are in excellent agreement, with the only discrepancies apparent in the twelfth decimal place or better; note that even at machine precision, the last several digits are often still fluctuating. A comparison of the convergence rates for the systems given by Eqs. 6 and 9 is shown in Fig. 2. The systems converge in a similar fashion. V. Large Scale Test Cases Several large scale test cases are used to evaluate the impact of the new scheme. For the cases shown, relative mea- sures in performance are used to establish benefits of the current approach. The meshes used here are not necessarily optimal for resolving the flowfields; the intent is solely to demonstrate algorithmic performance on a given discrete model. For the same reason, vehicle performance benefits obtained here through optimization are not claimed to be representative of what may be achieved for grid-resolved computations. Each of the grids shown has been obtained using the methods of Refs. 42 and 43. Note that although the cases described here involve compressible flows, an in- compressible implementation29,44 based on the method of artificial compressibility45 is also maintained. To show the impact of the current formulation, Table 3 lists timings for the major software components necessary for aerodynamic shape optimization. Operations associated with the surface grid evaluation and sensitivities are omit- ted; these costs are negligible in comparison with the other components. The computations for each configuration have been performed on different hardware with a range of compilers; timings for a given test case should only be compared relative to each other and not to other test cases. For each of the computations, the residuals for Eqs. 5 and Table 3 Timings for analysis and sensitivity analysis components. Wallclock Time, mins Test Case Analysis Sensitivity Analysis (Number of Mesh Adjoint Design Variables) Mesh Flowfield Flowfield Solution and Movement Solution Adjoint Solution Sensitivities Inviscid Stowed Morphing Vehicle 4 18 42 4 (175) Turbulent Cruise Morphing Vehicle 3 45 70 4 (150) Turbulent Wing–Body 8 228 184 10 (172) 5 American Institute of Aeronautics and Astronautics 9 have been reduced by ten orders of magnitude. All of the values shown in the table use an objective function based on the drag coefficient. The data shown in the table for each test case is discussed in more detail below. The time required for a flowfield adjoint solution is discussed at length in Ref. 31. The relative cost of these com- putations can vary widely depending on the coupling used for the turbulence closure equation. The storage strategy used for the exact linearization of the residual required for the adjoint solution can also have a large impact on the ef- ficiency of the adjoint solver. Here, the nearest neighbor terms are stored, while the larger stencil contributions are re- computed at each time step in order to save memory. If sufficient memory is available, the entire linearization may be stored and the relative cost of the adjoint solution for the flowfield can be drastically reduced by computing an ex- plicit matrix–vector product for the residual of Eq. 3. In order to save computational time during optimizations, the flow solver and its associated adjoint solver are typi- cally restarted from previous solutions, so that, in general, these solutions require less time as the optimization con- verges. Therefore the data indicated in Table 3 cannot necessarily be used to directly infer the cost of an entire opti- mization. However, the mesh linearizations are recomputed for every sensitivity analysis so the savings of the current adjoint formulation over the direct approach scale directly with the number of sensitivity analyses required over the course of an optimization. This benefit will be demonstrated in the following sections. Morphing Aircraft Configuration The first two test cases are based on an unmanned combat air vehicle being pursued by Lockheed-Martin and DARPA and is described in Ref. 46. The aircraft is designed to cruise in a conventional configuration and to fold its wings into a stowed position for the dash/attack portion of its mission. The first test case is based on the stowed configuration and is computed using the Euler equations. The grid used for this case is shown in Fig. 3, and contains 430,732 nodes and 2,543,772 tetrahedra. The freestream Mach number is 0.80 and the angle of attack is 2 degrees. For this test, the geometry has been parameterized using the package out- lined in Ref. 38 to create a set of 175 design variables describing the wing and fuselage. Timing results obtained on 16 Pentium IV processors for the various components needed for optimization are shown in Table 3. Here, the flow solver and its adjoint counterpart have been run 300 time steps. The analysis por- tion of the computation, consisting of a mesh movement and a flowfield solution, requires approximately 22 minutes of wallclock time. With a single solution of Eq. 9, the sensitivity analysis for all 175 design variables takes 46 min- utes. Note that the previous method based on direct differentiation would require 175 solutions to Eq. 6, or nearly 12 hours more wallclock time. This savings is the equivalent of 39 additional flowfield solutions. To demonstrate an optimization for the current test case, an objective function based on drag has been used in con- junction with an explicit lift constraint. In this manner, separate adjoint solutions are required to evaluate the sensitiv- ities of both the objective and the constraint. A subset of 63 of the 175 total shape variables are free to change, subject to side constraints that prohibit thinning of the geometry. The optimization package described in Ref. 47 is used to drive the design problem. The results of the optimization are shown in Fig. 4, where the drag coefficient has been re- duced by 52% and the lift coefficient satisfies the requested constraint value. For this case, the optimizer requested 26 Lift Drag LiftConstraint nt ent Coefficie D CD = –52% Coeffici Lift Drag 5 10 15 20 25 FunctionEvaluation Figure 3 Surface grid for morphing aircraft in dash Figure 4 Lift and drag coefficients during configuration. constrained optimization of dash configuration. 6 American Institute of Aeronautics and Astronautics o Rati D L⁄D = +24% D L/ 0 5 10 15 20 25 30 35 FunctionEvaluation Figure 5 Surface grid for morphing aircraft in cruise Figure 6 Lift-to-drag ratio during unconstrained configuration. optimization of cruise configuration. sensitivity analyses during the course of the optimization. Based on the timings above, the current adjoint formulation results in a savings of roughly 109 hours of wallclock time — or the equivalent of 364 flowfield solutions — in com- puting the mesh linearizations for the 63 variables. The next case is a turbulent flow computation for the baseline cruise configuration as shown in Fig. 5. The grid contains 534,525 nodes and 3,159,677 tetrahedra. Here, the freestream Mach number is 0.80, the angle of attack is 2 degrees, and the Reynolds number is 5 million based on the mean aerodynamic chord. This case has been run on 64 Opteron processors, and results for the component timings are again shown in Table 3. For these computations, the flow solver and associated adjoint solver have been run for 800 time steps using a loosely coupled formulation for the turbulence equation. The mesh linearizations for all 150 design variables have been computed in approximately four minutes. With the direct approach, equivalent results would require 150 solutions to Eq. 6 at an expense of roughly 450 minutes or the equivalent of 10 flowfield solutions. For this test case, the optimization goal is to perform an unconstrained maximization of the lift-to-drag ratio (L/D). The optimization package used for this computation is described in Ref. 48. Results for the L/D maximization are shown in Fig. 6, where a subset of 54 of the 150 shape design variables are allowed to change, subject to side con- straints that the geometry does not decrease in thickness. The final value of L/D is approximately 24% higher than its initial value. Here, the optimizer required seven sensitivity analyses. The present adjoint implementation for the mesh-related terms has yielded an estimated savings of approximately 19 hours of wallclock time, or the cost of 25 flowfield solutions, over the direct approach. Slotted Wing–Body Configuration The final example is fully turbulent flow over a transport wing–body configuration with a wing incorporating an advanced slotted airfoil concept. The mesh for this case is shown in Fig. 7 and contains 1,326,150 nodes and 7,744,304 tetrahedra. The freestream Mach number is 0.87, the angle of attack is 2.7 degrees, and the Reynolds num- ber is 3 million based on the mean aerodynamic chord. The computations shown in Table 3 have been performed using 64 Opteron processors. Here, the flow solver and its adjoint counterpart have each been run 1000 time steps, using the tightly coupled implementation of Ref. 31 for the turbulence equation. For this reason, the relative expense of the flowfield adjoint solution is less than that of the previous examples, since the complete Jacobian of the residual operator is required at every time step during the flow solution. The surface grids for this computation have been pa- rameterized using the method of Ref. 37 to create a set of 172 design variables describing the main wing and flap ge- ometries. For this case, the analysis problem takes 236 minutes of wallclock time versus 194 minutes for the sensitiv- ity analysis. The present adjoint formulation for the mesh terms yields a savings of approximately 23 hours in computing the sensitivities for all 172 variables. This is the equivalent of six additional flowfield solutions. An unconstrained maximization of L⁄D is performed for the current configuration, using a subset of 41 of the shape parameters in addition to the angle of attack. As in the previous examples, the geometry is constrained to main- tain its original thickness. The net result of the optimization procedure is a 32% increase in L/D as shown in Fig. 8. 7 American Institute of Aeronautics and Astronautics D L⁄D = +32% o ati R D L/ 0 5 10 15 20 25 30 35 FunctionEvaluation Figure 7 Surface grid for slotted wing–body Figure 8 Lift-to-drag ratio during unconstrained configuration. optimization of wing–body configuration. For this case, the optimizer required 26 sensitivity analyses during the course of the computation. The previous for- mulation for the grid sensitivities would have required an estimated additional 142 hours of wallclock time to per- form the same optimization. This savings is the equivalent of 37 flowfield solutions. VI. Summary and Conclusions An adjoint formulation to efficiently account for mesh sensitivities has been developed and implemented. In com- bination with an adjoint procedure for the aerodynamic flowfield, the new approach effectively removes the depen- dence on the size of the design space for a rigorous discrete sensitivity analysis. Unlike previous methods that re- quired explicit computation of grid sensitivity terms for each design variable, the current approach accounts for these terms through the solution of a single adjoint problem, equivalent in cost to the mesh movement scheme being used. The method has been implemented within a three-dimensional unstructured grid framework, and the resulting sensi- tivity derivatives have been shown to be in excellent agreement with a direct mode of differentiation using complex variables, as well as a previous implementation using handcoded direct differentiation for the mesh terms. Timings for the major components have been shown for several large-scale configurations, with dramatic savings obtained through the use of the adjoint formulation for the mesh linearizations. Demonstration optimizations yield improve- ments in vehicle performance for the discrete models used. Acknowledgments The authors wish to thank Bill Jones for providing the meshes used in this work and Jamshid Samareh for his assis- tance in parameterizing the various configurations. Results from Clyde Gumbert for the stowed morphing configura- tion are also appreciated. References 1Hicks, R.M., and Henne, P.A., "Wing Design by Numerical Optimization,"Journal of Aircraft, Vol. 15, 1978, pp. 407-412. 2Joh, C.-Y., Grossman, B., and Haftka, R.T., "Design Optimization of Transonic Airfoils,"Engineering Optimization, Vol. 21, 1993, pp. 1-20. 3Vanderplaats, G.N., Hicks, R.N., and Murman, E.M., "Application of Numerical Optimization Techniques to Airfoil Design," NASA Conference on Aerodynamic Analysis Requiring Advanced Computers," NASA SP-347, Part II, March 1975. 4Baysal, O., and Eleshaky, M.E., “Aerodynamic Sensitivity Analysis Methods for the Compressible Euler Equations,”Journal of Fluids Engineering, Vol. 113, 1991, pp. 681-688. 5Borggaard, J.T., Burns, J., Cliff, E.M., and Gunzburger, M.D., “Sensitivity Calculations for a 2-D Inviscid Supersonic Fore- body Problem,”Identification and Control Systems Governed by Partial Differential Equations, SIAM Publications, Philadelphia, 1993, pp. 14-24. 8 American Institute of Aeronautics and Astronautics 6Burgreen, G.W., and Baysal, O., “Aerodynamic Shape Optimization Using Preconditioned Conjugate Gradient Methods,” AIAA Paper 93-3322, 1993. 7Hou, G. J.-W., Maroju, V., Taylor, A.C., and Korivi, V.M., “Transonic Turbulent Airfoil Design Optimization with Automatic Differentiation in Incremental Iterative Forms,” AIAA 95-1692, 1995. 8Newman, J.C., and Taylor, A.C., “Three-Dimensional Aerodynamic Shape Sensitivity Analysis and Design Optimization Using the Euler Equations on Unstructured Grids,” AIAA 96-2464, 1996. 9Sherman, L.L., Taylor, A.C., Green, L.L., Newman, P.A., Hou, G.J.-W., and Korivi, V.M., “First- and Second-Order Aerody- namic Sensitivity Derivatives via Automatic Differentiation with Incremental Iterative Methods,” AIAA 94-4262, 1994. 10Young, D.P., Huffman, W.P., Melvin, R.G., Bieterman, M.B., Hilmes, C.L., and Johnson, F.T., “Inexactness and Global Con- vergence in Design Optimization,” AIAA 94-4386, 1994. 11Anderson, W.K., Newman, J.C., Whitfield, D.L., and Nielsen, E.J., "Sensitivity Analysis for the Navier-Stokes Equations on Unstructured Meshes Using Complex Variables,"AIAA Journal, Vol. 39, No. 1, 2001, pp. 56-63. 12Lyness, J.N., “Numerical Algorithms Based on the Theory of Complex Variables,” Proc. ACM 22nd Nat. Conf., Thomas Book Co., Washington, D.C., 1967, pp. 124-134. 13Lyness, J.N. and Moler, C.B., “Numerical Differentiation of Analytic Functions,”SIAM Journal of Numerical Analysis, Vol. 4, 1967, pp. 202-210. 14Martins, J.R.R.A., Kroo, I.M., and Alonso, J.J., "An Automated Method for Sensitivity Analysis using Complex Variables," AIAA 2000-0689, 2000. 15Newman, J.C., Anderson, W.K., and Whitfield, D.L., “Multidisciplinary Sensitivity Derivatives Using Complex Variables,” Mississippi State University Report No. MSSU-COE-ERC-98-08, 1998. 16Squire, W. and Trapp, G., “Using Complex Variables to Estimate Derivatives of Real Functions,”SIAM Review, Vol. 10, No. 1, 1998, pp. 110-112. 17Anderson, W.K. and Bonhaus, D.L., "Airfoil Design on Unstructured Grids for Turbulent Flows,"AIAA Journal, Vol. 37, No. 2, 1999, pp. 185–191. 18Anderson, W.K. and Venkatakrishnan, V., "Aerodynamic Design Optimization on Unstructured Grids with a Continuous Ad- joint Formulation,"Computers and Fluids, Vol. 28, No. 4, 1999, pp. 443–480. 19Elliott, J., "Discrete Adjoint Analysis and Optimization with Overset Grid Modelling of the Compressible High-Re Navier- Stokes Equations," 6th Overset Grid and Solution Technology Symposium, Fort Walton Beach, FL, Oct. 2002. 20Giles, M.B., Duta, M.C., Muller, J.-D., and Pierce, N.A., "Algorithm Developments for Discrete Adjoint Methods," AIAA Journal, Vol. 41, No. 2, February 2003, pp. 198-205. 21Iollo, A., Salas, M.D., and Ta’asan, S., "Shape Optimization Governed by the Euler Equations Using an Adjoint Method," ICASE Report No. 93-78, November 1993. 22Jameson, A., Pierce, N.A., and Martinelli, L., "Optimum Aerodynamic Design Using the Navier-Stokes Equations," AIAA 97- 0101, January 1997. 23Kim, C.S., Kim, C., and Rho, O.H., "Sensitivity Analysis for the Navier-Stokes Equations with Two-Equation Turbulence Models,"AIAA Journal, Vol. 39, No. 5, May 2001, pp. 838–845. 24Kim, H.-J., Sasaki, D., Obayashi, S., and Nakahashi, K., "Aerodynamic Optimization of Supersonic Transport Wing Using Unstructured Adjoint Method,"AIAA Journal, Vol. 39, No. 6, June 2001, pp. 1011–1020. 25Mohammadi, B., "Optimal Shape Design, Reverse Mode of Automatic Differentiation and Turbulence," AIAA 97-0099, Janu- ary 1997. 26Nemec, M. and Zingg, D.W., "Towards Efficient Aerodynamic Shape Optimization Based on the Navier-Stokes Equations," AIAA Paper 2001-2532, 2001. 27Newman III, J.C., Taylor III, A.C., and Burgreen, G.W., "An Unstructured Grid Approach to Sensitivity Analysis and Shape Optimization Using the Euler Equations," AIAA 95-1646, 1995. 28Nielsen, E.J., "Aerodynamic Design Sensitivities on an Unstructured Mesh Using the Navier-Stokes Equations and a Discrete Adjoint Formulation," Ph.D. Dissertation, Dept. of Aerospace and Ocean Engineering, Virginia Polytechnic Inst. and State Univ., December 1998. 29Nielsen, E.J. and Anderson, W.K., "Aerodynamic Design Optimization on Unstructured Meshes Using the Navier-Stokes Equations,"AIAA Journal, Vol. 37, No. 11, 1999, pp. 1411–1419. 30Nielsen, E.J. and Anderson, W.K., "Recent Improvements in Aerodynamic Design Optimization on Unstructured Meshes," AIAA Journal, Vol. 40, No. 6, 2002, pp. 1155–1163. 31Nielsen, E.J., Lu, J., Park, M.A., and Darmofal, D.L., "An Implicit, Exact Dual Adjoint Solution Method for Turbulent Flows on Unstructured Grids,"Computers and Fluids, Vol. 33, No. 9, 2004, pp. 1131-1155. 32Reuther, J.J., Jameson, A., Alonso, J.J., Rimlinger, M.J., and Saunders, D., "Constrained Multipoint Aerodynamic Shape Opti- mization Using an Adjoint Formulation and Parallel Computers,"Journal of Aircraft, Vol. 36, No. 1, 1999, pp. 51–60. 33Soemarwoto, B., "Multipoint Aerodynamic Design by Optimization," Ph.D. Dissertation, Dept. of Theoretical Aerodynamics, Delft University of Technology, December 1996. 34Soto, O. and Lohner, R., "A Mixed Adjoint Formulation for Incompressible Turbulent Problems," AIAA 2002-0451, 2002. 35Sung, C. and Kwon, J.H., "Aerodynamic Design Optimization Using the Navier-Stokes and Adjoint Equations," AIAA 2001- 0266, 2001. 36Sundaram, P., Agrawal, S., and Hager, J.O., "Aerospace Vehicle MDO Shape Optimization Using ADIFOR 3.0 Gradients," AIAA 2000-4733, 2000. 37Samareh, J.A., "A Novel Shape Parameterization Approach," NASA TM-1999-209116, May 1999. 38Samareh, J.A., "Aerodynamic Shape Optimization Based on Free-Form Deformation," AIAA 2004-4630, 2004. 39Saad, Y. and Schultz, M.H., "GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Sys- tems,"SIAM Journal of Scientific and Statistical Computing, Vol. 7, No. 3, 1986, pp. 856–869. 9 American Institute of Aeronautics and Astronautics 40Schmitt, V. and Charpin, F., "Pressure Distributions on the ONERA M6 Wing at Transonic Mach Numbers," Experimental Database for Computer Program Assessment, AGARD-AR-138, May 1979, pp. B1-1-B1-44. 41Kleb, W.L., Nielsen, E.J., Gnoffo, P.A., Park, M.A., and Wood, W.A., "Collaborative Software Development in Support of Fast Adaptive AeroSpace Tools (FAAST)," AIAA 2003-3978, 2003. 42Jones, W.T., "GridEx - An Integrated Grid Generation Package for CFD," AIAA 2003-4129, 2003. 43Marcum, D.L., and Weatherill, N.P., "Unstructured Grid Generation Using Iterative Point Insertion and Local Reconnection," AIAA Journal, Vol. 33, No. 9, 1995, pp. 1619-1625. 44Anderson, W.K., Rausch, R.D., and Bonhaus, D.L., "Implicit/Multigrid Algorithms for Incompressible Turbulent Flows on Unstructured Grids,"Journal of Computational Physics, Vol. 128, 1996, pp. 391-408. 45Chorin, A.J., "A Numerical Method for Solving Incompressible Viscous Flow Problems,"Journal of Computational Physics, Vol. 2, 1967, pp. 12-26. 46Love, M.H., Zink, P.S., Stroud, R.L., Bye, D.R., and Chase, C., "Impact of Actuation Concepts on Morphing Aircraft Struc- tures," AIAA 2004-1724, 2004. 47Gill, P.E., Murray, W., Saunders, M.A., and Wright, M.H., "User’s Guide for NPSOL 5.0: A FORTRAN Package for Nonlin- ear Programming," Technical Report SOL 94, 1995. 48Kaufman, L. and Gay, D., "PORT Library: Optimization and Mathematical Programming - User’s Manual," Bell Laboratories, 1997. 10 American Institute of Aeronautics and Astronautics