Development of an efficient aerodynamic shape optimization framework

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 79 (2009) 2373–2384

Development of an efficient aerodynamic shape optimization framework Jong-Eun Kim ∗ , Vinay N. Rao, Roy P. Koomullil, Doug H. Ross, Bharat K. Soni, Alan M. Shih Department of Mechanical Engineering, University of Alabama at Birmingham, Birmingham, AL 35294-4461, USA Available online 31 January 2009

Abstract Although many efforts have been made to develop an aerodynamic shape optimization (ASO) framework, iterative grid generation of the complex configuration within the optimization loop has still been a critical barrier. In this paper, an efficient ASO framework is developed by integrating a parametric grid generator, an optimization toolkit, and a flow solver. A geometry-grid template toolkit is developed to address the need to produce a large number of grids in a timely manner for the parametric design study. An objectoriented optimization toolkit that allows a flexible and extensible interfacing with user-specific codes is used. An in-house full Navier–Stokes flow solver is developed and used in the framework. Code integration is achieved using a black-box interface with script files. Two ASO applications and their optimum solutions are presented to demonstrate the success of this framework. © 2009 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Aerodynamic shape optimization; Parametric grid generator; Optimization toolkit; Flow analysis; Surrogate-based optimization

1. Introduction With the help of high performance computing platforms and efficient numerical algorithms for flow analysis and optimization, aerodynamic shape optimization (ASO) techniques have been applied to improve the aerodynamic performances of aerospace vehicles. A survey paper by Newman et al. [13] provides a comprehensive discussion of aerodynamic shape and design sensitivity analysis and optimization, based on advanced computational fluid dynamics (CFD). In ASO, research has focused on sensitivity analysis and shape parameterization techniques. For efficient sensitivity analysis, continuous or discrete adjoint formulations have been used, in which the necessary gradients are obtained through the solution of an adjoint system of the governing equations of interest [7,14]. For multidisciplinary optimization of aerospace applications, a variety of shape parameterization techniques have been developed such as basic vector, domain element, partial differential equation, discrete, analytical, polynomial and spline, free form deformation, and computer aided design (CAD)-based approaches [17]. Among them, CAD-based parameterization [12] has been popularly used for ASO because of its capability of creating details of three-dimensional complex configurations. Although this approach has made considerable progress, its practical application is still hindered by significant computational ∗

Corresponding author. BEC 257, 1530 3rd Ave. S., Birmingham, AL 35294-4461, USA. Tel.: +1 205 975 5889; fax: +1 205 975 7244. E-mail addresses: [email protected] (J.-E. Kim), [email protected] (V.N. Rao), [email protected] (R.P. Koomullil), [email protected] (D.H. Ross), [email protected] (B.K. Soni), [email protected] (A.M. Shih). 0378-4754/$36.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2009.01.012

2374

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

time for iterative grid generation that must reflect the change of configuration within the optimization loop. Moreover, most CAD software packages do not offer grid generation tools, and most grid generators do not have the capability to create sophisticated geometric configurations that CAD packages offer. In this study, a CAD-based parametric grid generator called “MiniCAD” is developed to address the need to produce a large number of grids in a timely manner for the parametric design study. Functionalities for geometry creation and modification are defined using NURBS (non-uniform, rational B-spline) curve and surface and volume formulations. An optimization toolkit, Design and Analysis Kit for Optimization and Terascale Applications (DAKOTA), is used as an optimization framework. This toolkit provides a flexible and extensible interfacing with user-specific codes. Instead of directly operating a computationally expensive flow model in the optimization loop, surrogate-based optimization methodology is applied for the ASO applications. A flow simulation system that is used for the current study is developed for a generalized grid framework. Code integration is achieved by a so-called black-box interface using script files in which the user does not need to have access to the source codes. 2. Development of an efficient ASO framework 2.1. Optimization toolkit Optimization is a mathematical process to obtain an optimal design variable set that achieves design objectives by minimizing or maximizing objective functions. In addition, all design constraints must be satisfied. In ASO, the drag coefficient or the ratio of lift and drag can be defined as an objective function. Shape parameters that affect aerodynamic performances at large can be defined as design variables. An optimization toolkit named DAKOTA, developed at Sandia National Laboratories [2], is used for this study. It was originally developed as a collection of gradient-based and non-gradient-based optimization software and now includes methods for statistical data sampling, surrogate modeling, uncertainty quantification, and a variety of other statistical and mathematical software tools useful to design engineers. DAKOTA has a generic interface capability that permits to couple it to a variety of commercial and custom-developed engineering codes. 2.2. Parametric grid generator A geometry-grid generation template tool, MiniCAD [18], has been developed to address the need to produce a large number of grids in a timely manner for parametric design study. It provides relatively sophisticated CAD capabilities that are generally not available in a grid generation tool. It also allows parametric geometry generation, which is a critical function in this design optimization framework. Once a template is defined, a user can modify values for multiple geometric parameters either interactively or through scripting using standard Python language. The new geometry and its corresponding mesh will be produced automatically. When the configuration changes, due to the design need, the template can be quickly modified to reflect the change. Functionalities for geometry creation and modification are defined using NURBS curve, surface, and volume formulations [19]. A NURBS curve of order k is defined as n 

C(u) =

wi di Nik (u)

i=0 n 

(1) wi Nik

(u)

i=0

where di (i = 1, . . . , n) denotes the deBoor control polygon and wi are the non-equal weights associated with each control points. When wi = 1.0, Eq. (1) represents a B-Spline curve. The Nik (u) is the normalized B-Spline basis function of order k and is defined as Nik (u) =

k−1 (u) (u − Ti ) Nik−1 (u) (Ti+k − u) Ni+1 + Ti+k−1 − Ti Ti+k − Ti+1

(2)

Ni1 (u) =

if Ti ≤ u ≤ Ti+1 0, otherwise

(3)

1,

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

2375

From the formulation of a NURBS curve, a tensor product NURBS surface with order k1 and k2 is shown as m  n 

S(s, t) =

wij dij Nik1 (s) Njk2 (t)

i=0 j=0 m  n 

(4) wij Nik1 (s) Njk2 (t)

i=0 j=0

where dij (i = 0, . . . , m, j = 0, . . . , n) are 3D NURBS control net and wij are the weights associated with each control point. Nik1 (u) and Njk2 (u) denote the normalized basis functions of order k1 and k2 over two knot vector T1 = Ti (i = 0, . . . , m + k1 ) and T2 = Tj (j = 0, . . . , n + k2 ), respectively. The basis function is defined recursively the same as the one defined in the NURBS curve. Taking the parametric forms in NURBS formulations, one can modify the geometry accurately. Therefore, in this study, NURBS-based MiniCAD is used to provide the above geometry-grid generation template capabilities with Python scripts to serve as the simulation process control. Fig. 1 shows a parametric duct geometry with different shapes of cross-section. The updated design variables calculated from the optimization toolkit are fed back to the MiniCAD template to regenerate the geometry and the associated meshes automatically. Five software components are developed and integrated to form the MiniCAD framework. (1) User interaction environment: the user interaction environment is a set of Python-based software modules that forms the Praxis Environment. The Praxis Environment provides a consistent set of tools for application development and inter-application connectivity. (2) Geometry design and CAD data interface: geometry-related functionalities are mainly maintained in the GeometryGrid Toolkit (GGTK) library [3]. GGTK is a standalone software library of functions for geometric modeling, topological representation, and grid generation. It serves as the underlying library for the MiniCAD template framework. (3) Grid generation algorithms: grid generation using structured [20], unstructured [5], or hybrid [6] grid topology is available in MiniCAD framework. (4) Grid pattern library: grid patterns are topological structures that include grid generation information. They can be reused by connection to different geometries. This feature is unique to MiniCAD. The benefit of grid patterns is that a user can map a grid pattern onto a topologically equivalent geometry to quickly create the grid. (5) Template library: creation of template parameters includes the writing of equations to form the interrelationships between object attributes and the parameters. After the template is completed, parameters are modified and the update operation is initiated to execute the equations and the object’s parametric relationships to recreate the geometry and grid based on the new parameter values. 2.3. Flow simulation system A flow simulation system named HYB3D that is used for the current study is developed for a generalized grid framework, in which the discretization of the physical domain can be of structured, unstructured, or an agglomeration of cells with an arbitrary number of faces [9–11]. The integral form of the Navier–Stokes equations is taken as the governing equations for the fluid flow and can be written as    ∂Q d + F (Q) · n ds = F v (Q) · n ds (5)  ∂t ∂ ∂

where Q is the conserved variable vector, F (Q) is the inviscid flux vector, F v (Q) is the viscous flux vector,  is the control volume, ∂ is the control surface, n is the outward pointing unit normal to the control surface, and ds is the elemental area. A cell-centered finite volume scheme, in which the cell averaged values are assumed to be at the cell center, is used for the discretization of the governing equations given in Eq. (5). The discretized form of the

2376

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

Fig. 1. MiniCAD system for a duct geometry. (a) Triangular cross-section and (b) square cross-section.

Navier–Stokes equations can be written as k k   Q Vi = − Fij · n dsj + Fijv · n dsj t j=1

(6)

j=1

where Vi is the volume of the cell i, j represents the faces that form the cell, and Fij · n ds and Fijv · n ds are the inviscid and viscous numerical fluxes crossing the jth face of the cell i, respectively. The inviscid numerical fluxes passing through the cell faces are calculated by Roe’s approximate Riemann solver [15] as an exact solution for a linearized

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

2377

Riemann problem. Using the approximate Riemann solver, the flux through a cell face is given by Fij =

1 ¯ [F (QL ) + F (QR ) − |A|(Q R − QL )] 2

(7)

¯ is the Roe-averaged matrix [15]. where R and L represent cells on the left and right sides of the cell-face, and |A| The flow properties on either sides of the cell-face for the Riemann problem is taken as the cell averaged values for first order scheme, while a linear reconstruction of the flow variables is used for the second order scheme. The linear reconstruction of the flow variables is achieved using the Taylor’s series expansion and the flow properties in the neighboring cells. The gradient of the flow variables at the cell center, required for the Taylor’s series expansion, is estimated using Green’s theorem with the control volume taken as the cell itself or using the principles of least squares [9]. The values of the conserved variable vector at the nodes are estimated based on a weighted-average procedure. During the reconstruction process, local extrema may be created, which may produce spurious values in regions where there are sharp jumps in the flow variables such as shocks, contact discontinuities, expansion regions, etc. To avoid the generation of spurious values, a limiter function is employed as Q(x, y, z) = Q(xi , yi , zi ) + φ(∇Q)(xi ,yi ,zi ) · r + O(r 2 )

(8)

where ∇Q is the gradient of the vector Q, r is the vector pointing from the centroid of the cell to the centroid of the face, and φ is the limiter function and its value is limited between zero and one. Limiters by Venkatakrishnan [22], and Barth and Jesperson [1] are implemented in the generalized grid framework. In the case of implicit schemes, the numerical flux crossing the cell face is a function of the conserved variables at the (n + 1)th time level. The flux vector is linearized before the evaluation of the flux crossing the cell faces. The resulting linear system can be written as ⎤ ⎡ n  k  k    ∂H ∂Hij n V ij ⎦ Qni + ⎣ iI+ Qnj = −ni (9) t ∂Qi ∂Qj j=1

j=1

where Hij = Fij · n ds. The Jacobian matrices (∂Hij /∂Qi ) and (∂Hij /∂Qj ) can be estimated using approximate analyt¯ to be constant, or by a numerical approach [23]. The block sparse ical Jacobians, by taking the Roe-averaged matrix |A| matrix system resulting from the linearization of the governing equations is solved using a symmetric Gauss–Seidel algorithm. Newton iterations [23] are used to improve the temporal accuracy for the case of time-accurate simulations using the implicit scheme. The parallelization of the code is achieved by decomposing the physical domain into different regions using METIS [8], which utilizes the graph of the grid to perform the decomposition. Message passing interface (MPI) [4] is used to pass off the information across the block interfaces. 2.4. Code integration In DAKOTA, external grid generators and simulation codes can be coupled by a direct interface or a black-box interface. Although direct interface is more computationally efficient than black-box interface, it involves significant code development time. In black-box interface, DAKOTA and external codes remain entirely independent, with data transferred through writing and reading text files. This interface scheme makes it easy to build, control and utilize the parallel computing capabilities that simulation codes have. In this study, code integration is achieved using the black-box interface where the user does not need to have access to the source codes. The entire optimization process is automatically executed using script files written by PERL, BASH, and C for data extraction, manipulation, file handling, and code execution. Fig. 2 shows a flow diagram for the ASO framework. The DAKOTA input file contains a description of design variables (variable types, initial points, and lower and upper bounds), optimization algorithm and strategy, interface method, and responses (numbers of objective functions and constraints, and method for gradient calculation). This input file is executed on the command line of a “job” file to invoke the process on parallel processors. DAKOTA generates a new set of design variables based on selected optimization algorithm and strategy. This new set of design variables is parsed through the input file of the MiniCAD to produce a new grid. An in-house filter is required to convert the grid file format into an input file format of the specific CFD code. Flow responses of interest are extracted from the output files of the flow solver. The objective function is

2378

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

Fig. 2. Flow diagram for the ASO framework.

Fig. 3. Procedure of surrogate-based optimization.

computed with the flow responses and is fed to DAKOTA to continue the process. The entire process is automated and executed in parallel. The file parsing and data transfer are achieved by PERL and BASH scripts. Template input files are used to create new input files. Coarse-grained parallelism in DAKOTA is utilized to greatly minimize computational time [2]. 2.5. Optimization strategies Due to the requirement of many iterations to obtain an optimum solution, a general optimization scheme is not feasible for the CPU-intensive ASO. Instead of directly operating a computationally expensive flow model in the optimization loop, surrogate-based optimization (SBO) methodology can be applied. The SBO utilizes data sampling and surface fitting strategies to build a computationally efficient mathematical model for the purpose of numerical optimization. Fig. 3 illustrates the entire process of the SBO. In order to construct a surface fitting model, a data sampling process is necessary to select certain sets of design variables with which analyses will be performed using a simulation code. In this ASO, several sets of shape parameters are generated for MiniCAD to prepare the geometry and grids. The number of points is dependent upon the number

Fig. 4. LHS in a 2D design space.

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

2379

Fig. 5. Shape parameters for a generic rocket. Table 1 Design points and calculated objective function. No.

x1

x2

x3

x4

x5

Cd

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

20.857 28.891 27.833 37.207 35.655 31.052 38.485 22.409 35.847 25.489 25.484 22.654 29.211 39.530 35.658 32.197 36.968 21.050 28.469 27.528

0.119 0.098 0.107 0.101 0.129 0.127 0.108 0.102 0.096 0.119 0.099 0.127 0.106 0.095 0.114 0.098 0.100 0.125 0.119 0.114

158.653 160.587 163.599 152.847 155.772 158.870 163.105 160.038 150.272 153.521 157.340 154.212 152.690 163.277 161.551 156.065 163.040 154.239 152.023 160.777

135.120 131.437 138.701 133.083 137.648 139.139 135.113 136.222 130.260 133.674 136.160 138.515 132.011 134.069 130.257 135.304 137.416 139.522 133.766 130.900

0.047 0.038 0.040 0.043 0.032 0.049 0.034 0.040 0.045 0.031 0.039 0.037 0.048 0.033 0.043 0.041 0.047 0.032 0.044 0.037

0.0124712 0.0266225 0.0233250 0.0292298 0.0298239 0.0253959 0.0376303 0.0125353 0.0227407 0.0137885 0.0129446 0.0118714 0.0156512 0.0384298 0.0462965 0.0223127 0.0326135 0.0096642 0.0159562 0.0274988

of design variables and the types of data sampling methods. DAKOTA provides many sampling methods such as Latin hypercube, Monte Carlo, orthogonal array, central composite and factorial samplings. Among them, the Latin hypercube sampling (LHS) is an attractive method when a large number of design variables are needed. This method allows the user to determine the number of samples. Fig. 4 shows a LHS sampling for a 2D design space in a case of having four design points. The stars are sample sites randomly selected inside each bin. Surface fitting is the process of constructing global approximation to system behavior using a proper polynomial or statistical functions. This is based on the simulation results calculated at selected design points through the data sampling process. This approach is useful when the functional relationship between some responses (e.g. drag coefficient) and independent input variables (e.g. shape parameters) is extremely complicated. Polynomial regression, the most popular surface fitting method, is to use a proper polynomial [16]. Consider a response variable, y, dependent on a number of variables, x. y = η(x)

(10)

Table 2 Initial design point and optimum designs.

Inital design Optimum (Kriging) Optimum (ANN)

x1

x2

x3

x4

x5

Cd

30.000 20.000 20.000

0.110 0.090 0.090

157.500 156.197 150.000

135.000 135.790 135.975

0.040 0.030 0.050

0.0232183 0.0044776 0.0049467

2380

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

The exact functional relationship can be approximated as η(x) ≈ f (x)

(11)

The approximating function, f, can be assumed to be a summation of basic functions f (x) =

L 

ai i (x)

(12)

i=1

Fig. 6. Optimal aerodynamic shape of a rocket. (a) Initial design, (b) optimal design (Kriging), and (c) optimal design (ANN).

Fig. 7. Shape parameters for a bent feedline pipe.

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

2381

where L is the number of basic functions, . The popular second order quadratic function approximation is expressed as n n n n     f (x) = a0 + ak xk + akk xk2 + akm xk xm (13) k=1

k=1

k=1 m=1,k
where n is the number of design variables. The constants, a, have to be determined in order to minimize the sum of the square error. ⎧⎡ ⎤2 ⎫ ⎪ P P ⎪ L ⎨ ⎬    ⎣yi (x) − (14) {[yi (x) − fi (x)]2 } = aj j (x)⎦ ⎪ ⎪ ⎩ ⎭ i=1 i=1 j=1 where P is the number of the design points and y is the exact functional response at the design point. Using the method of least squares and matrix algebra, the solution to the unknown coefficients is ¯ T X) ¯ −1 X ¯T y a = (X

(15)

DAKOTA provides other surface fitting methods such as Kriging interpolation, artificial neural network, and multivariate adaptive regression splines [2]. Program CONMIN [21], which is based on the method of feasible direction, is selected as an optimizer in this study. 3. Applications 3.1. Test case I: a generic rocket As a test case, the shape of a generic rocket is optimized to minimize the coefficient of drag force (Cd ) acting on the rocket surface. Five shape parameters (design variables) shown in Fig. 5 are used for the optimization purpose. The ranges of the design variables are set to 20 ≤ x1 ≤ 40, 0.09 ≤ x2 ≤ 0.13, 150 ≤ x3 ≤ 165, 130 ≤ x4 ≤ 140 and 0.03 ≤ x5 ≤ 0.05. The flow condition is taken as an inviscid flow with a free-stream Mach number of 0.6 and zero degree angle of attack. The grid used for this simulation consists of 373,626 elements, 403,534 nodes, and 1,149,942 faces. Twenty design points are selected through the LHS and listed in Table 1 with calculated Cd values. The optimum

Fig. 8. An optimal design of the feedline pipe.

2382

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

results through two different surface fitting methods, Kriging and artificial neural network (ANN), are shown in Table 2. The initial design point is set to the middle points of the range of each design variable. Fig. 6 illustrates the profile change from an initial design to optimal designs. The optimal designs using different surface fitting methods show a similar trend. From the results, the angle (x1 ) and the height (x2 ) that reach their lower bounds are most important design parameters to reduce the coefficient of drag force. The height (x5 ) has a minor effect. The optimum shape parameters significantly reduce the Cd value. 3.2. Test case II: a bent feedline pipe As another test case, the shape of a bent feedline pipe that has non-circular cross section is optimized. The objective is to maximize uniformity (by minimizing standard deviation, S.D.) of outlet flow velocity with five shape parameters shown in Fig. 7. The ranges of the design variables are set to 1.5 ≤ x1 ≤ 2.5, 0.5 ≤ x2 ≤ 1.5, 1.5 ≤ x3 ≤ 2.5, 0.5 ≤ x4 ≤ 1.5 and 2.0 ≤ x5 ≤ 4.0. The shape of the inlet and outlet is invariable and is an ellipse with major radius 2.0 and minor radius

Fig. 9. Velocity profile of outlet flow to check uniformity. (a) Initial design (S.D. = 13.68) and (b) optimal design (S.D. = 13.55).

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

2383

1.0. The inlet flow velocity is 68 m/s, and Reynolds number is set to 1000. The grid used for this simulation consists of 370,070 elements, 171,312 nodes, and 899,618 faces. Twenty-five design points are selected through the LHS. The quadratic polynomial regression is used for the surface fitting. Fig. 8 illustrates a design with optimal shape parameters (x1 = 2.5, x2 = 1.1, x3 = 1.5, x4 = 1.5 and x5 = 4.0). Fig. 9 shows the comparison of velocity uniformity of outlet flow between an initial design (with middle points of the ranges of each design variables) and the optimal design. The optimal design slightly reduces the value of S.D. from 13.68 to 13.55. 4. Conclusion An efficient aerodynamic shape optimization framework, which consists of an optimization toolkit, a parametric grid generator, a flow solver, and interfacing codes and scripts, was successfully developed and applied to some design problems. After completing the code integration, total time taken for each test case was approximately 3 days: 2 days for MiniCAD system for specific geometry and grid preparation, and less than 1 day for the execution of optimization including post-processing using multiprocessors. This framework could be applicable and extendable to any multidisciplinary optimization by coupling with any simulation code using box-based interfacing. Acknowledgement This research effort is supported in part by funding from the NASA Constellation University Institutes Project (CUIP) No. NCC3-994, and Department of Defense Programming Environment and Training (PET) program, CFDK4-02. References [1] T.J. Barth, D.C. Jespersen, The design and application of upwind schemes on unstructured meshes, AIAA Paper 89-0366, 1989. [2] M.S. Eldred, A.A. Giunta, B.G. van Bloemen Waanders, S.F. Wojtkiewicz, Jr., W.E. Hart, M.P. Alleva, DAKOTA Reference Manual: Version 3.2, Technical Report, Sandia National Laboratories, Albuquerque, NM, 2004 (online version available: http://endo.sandia.gov/DAKOTA). [3] S. Gopalsamy, D.H. Ross, Y. Ito, A.M. Shih, Structured Grid Generation over NURBS and Facetted Surface Patches by Reparameterization, in: Proceedings of the 14th International Meshing Roundtable, San Diego, CA, September 2005. [4] W. Gropp, E. Lusk, A. Skjellum, Using MPI, Portable Parallel Programming with the Message-Passing Interface, The MIT Press, 1999. [5] Y. Ito, P.C. Shum, A.M. Shih, B.K. Soni, K. Nakahashi, Robust generation of high-quality unstructured meshes on realistic biomedical geometry, International Journal for Numerical Methods in Engineering 65 (6) (2006) 943–973. [6] Y. Ito, A.M. Shih, B.K. Soni, K. Nakahashi, An approach to generate high quality unstructured hybrid meshes, AIAA Paper 2006-0530, 2006. [7] A. Jameson, L. Martinelli, S. Cliff, S. Thomas, Aerodynamic shape optimization of transonic and supersonic aircraft configurations, AIAA Paper 2005-1013, 2005. [8] G. Karypis, V. Kumar, METIS: Unstructured Graph Partitioning and Sparse Matrix Reordering System, Version 2.0, Department of Computer Science, University of Minnesota, August 1995. [9] R.P. Koomullil, D.S. Thompson, B.K. Soni, Iced airfoil simulation using generalized grids, The Journal of Applied Numerical Mathematics 46 (3–4) (2003) 319–330. [10] R.P. Koomullil, B.K. Soni, Flow simulation using generalized static and dynamics grids, AIAA Journal 37 (12) (1999) 1551–1557. [11] R.P. Koomullil, Flow Simulation System for Generalized Static and Dynamic Grids, Ph.D. Dissertation, Department of Aerospace Engineering, Mississippi State University, MS, May 1997. [12] M. Nemec, M.J. Aftosmis, T.H. Pulliam, CAD-based aerodynamic design of complex configurations using a cartesian method, AIAA Paper 2004-0113, 2004. [13] J.C. Newman, A.C. Taylor, R.W. Barnwell, P.A. Newman, G.J.W. Hou III, Overview of sensitivity analysis and shape optimization for complex aerodynamic configuration, Journal of Aircraft 36 (1) (1999) 87–96. [14] E.J. Nielsen, W.K. Anderson, Recent improvements in aerodynamic design optimization on unstructured meshes, AIAA Journal 40 (6) (2002) 1155–1163. [15] P.L. Roe, Approximate Riemann solvers, parameter vector, and difference schemes, Journal of Computational Physics 43 (1981) 357–372. [16] W.J. Roux, N. Stander, R.T. Haftka, Response surface approximations for structural optimization, International Journal for Numerical Methods in Engineering 42 (1998) 517–534. [17] J.A. Samareh, Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization, AIAA Journal 39 (5) (2001) 877–884. [18] A.M. Shih, D.H. Ross, M. Dillavou, S. Gopalsamy, B.K. Soni, A geometry-grid generation template framework for propellant delivery system, AIAA Paper 2006-5045, 2006.

2384

J.-E. Kim et al. / Mathematics and Computers in Simulation 79 (2009) 2373–2384

[19] A.M. Shih, T.-Y. Yu, S. Gopalsamy, Y. Ito, B.K. Soni, Geometry and mesh generation for high fidelity computational simulations using non-uniform rational B-splines, Applied Numerical Mathematics 55 (3) (2005) 368–381. [20] A.M. Shih, T.-Y. Yu, B.K. Soni, Interactive grid generation and NURBS applications, Journal of Applied Mathematics and Computation 65 (1–3) (1994) 279–294. [21] G.N. Vanderplaats, CONMIN—A Fortran Program for Constrained Function Minimization, User’s Manual, NASA TM X-62282, August 1973. [22] V. Venkatakrishnan, On the accuracy of limiters and convergence to steady state solutions, AIAA Paper 93-0880, 1993. [23] D.L. Whitfield, L. Taylor, Discretized Newton-relaxation solution of high resolution flux-difference split schemes, AIAA Paper 91-1539, 1991.