An Improved Parallel Optimization Framework for Transonic Airfoil Design

Variable Fidelity Optimization (VFO) is used to obtain a minimum drag transonic airfoil, at constant lift, subject to thickness and pitching moment constraints using several low fidelity solvers. VFO has emerged as an attractive method of performing, both, high-speed and high-fidelity optimization. VFO uses computationally inexpensive low-fidelity models, complemented by a surrogate to account for the difference between the high-and low-fidelity models, to obtain the optimum of the function efficiently and accurately. The authors' original Variable Fidelity (VF) framework is modified for increased efficiency and accuracy by incorporating parallel evaluation and more constraints to the optimization problem. The method is found to be efficient and capable of finding the optimum that closely agrees with the results of high-fidelity optimization alone.


INTRODUCTION
Given the high cost of CFD and optimization, a prominent area of research today is to find ways to reduce the computational time while retaining the highfidelity of the analysis. In the area of aerodynamic optimization, the variable fidelity (also called multifidelity) method has quickly grown in popularity (Alexandrov et al., 2001;Keane, 2003;Gano et al., 2005;Forrester et al., 2006;Forrester et al., 2007;Nelson et al., 2007;Lehner et al., 2010).
Variable-fidelity and other model management methods have been developed to solve optimization problems that involve simulations with large computational expense (Gano et al., 2005;Huang et al., 2006). In many engineering design problems, differing levels of fidelity can model the system of interest. Higher-fidelity models typically incorporate more detailed physics and are computationally expensive to evaluate than lower fidelity models. The lower-fidelity models are typically much cheaper to evaluate, but designs produced by using these models neglect important physical effects included in the more expensive higher-fidelity models. In aircraft design, Navier-Stokes and Euler equations are examples of two computational models with different fidelity.
Variable-fidelity optimization has emerged as an attractive method of performing, both, high-speed and high-fidelity optimization (Alexandrov et al., 2001;Keane, 2003;Leary et al., 2003;Gano et al., 2005;Forreste et al., 2006;Huang et al., 2006;Forrester et al., 2007;Nelson et al., 2007;Lehner et al., 2010). These algorithms attempt to leverage information from computationally inexpensive low-fidelity models to reduce the time required to converge to the optimum of the high-fidelity function. This is usually accomplished by building a computationally inexpensive surrogate for the high-fidelity model. The surrogate is, often, iteratively optimized and updated as new high-fidelity data become available. In variable-fidelity algorithms, the surrogate for the high-fidelity model usually consists of the low-fidelity model plus a correction term that models the difference between the high-and lowfidelity models, calibrated at selected sample points in the design space (Lehner et al., 2010).
The variable-fidelity method has its roots in the empirical model-building theory -the idea of endowing a surrogate with some discipline related properties to increase its accuracy-and the past three decades have seen rapid increase in its development and use (Simpson et al., 2008). Insightful reviews of the variable-fidelity method can be found in Huang et al. (2006) and Simpson et al. (2008). The method has improved immensely; from its initial form requiring Taylor polynomials (Alexandrov et al., 2001) to its current incarnation using a variety of modern surrogate models like Kriging and neural networks (Leary et al., 2003;Gano et al., 2005); and from using a variety of multiplicative, additive and hybrid corrections (Gano et al., 2005), to the method of Co-Kriging (Forrester et al., 2007). In the optimization context, the method has been used in both derivative-based (Alexandrov et al., 2001) and derivative free (Booker et al., 1999) optimization.
The authors have previously demonstrated a Variable-Fidelity Optimization (VFO) framework for the solution of a 2-dimensional aerodynamic design optimization problem (Zahir and Gao, 2011). The importance of prudently selecting the LF model is also highlighted: Inappropriate LF models mislead the optimization algorithm towards infeasible points in the design space (Zahir and Gao, 2011). The authors' framework worked well: The lift-to-drag ratio, K, was within 6% of the HF optimum and provided significant savings in computation time. Detailed results can be found in Zahir and Gao (2011). In the quest for greater efficiency and accuracy, the framework is modified to perform parallel computations. This study presents a minimum drag VFO of a transonic airfoil, at constant lift, subject to thickness and pitching moment constraints using several low fidelity solvers.
This study is arranged as follows: First the major technology pieces used in the design are described; next, the VFO framework is presented followed by the results of the optimization. The VF optimization results are also compared to the results of direct optimization, where the HF solver, alone, was coupled to the optimizer to find the optimum result.

MATERIALS AND METHODS
Flow solver: In this study, the RAE-2822 airfoil is chosen for optimization. An indigenously developed 2D compressible Navier-Stokes solver-using the LU-SGS time stepping scheme, the Roe upwind scheme and multigrid acceleration-capable of being used in, both, Euler and Reynolds-Averaged-Navier-Stokes (RANS) mode (Su et al., 2008) is used as the flow solver. Euler and low resolution (small grid) RANS solvers are used as the LF solvers, while RANS with the K-ω turbulence model is used as the HF solver.
The Euler and RANS solvers use a C-type mesh extending 20 chord lengths downstream of the trailing edge. The first grid line is 2x10 -6 units from the airfoil surface for the RANS solver and 0.01 units for the Euler solver. The HF solver uses a grid size of 216 cells around the airfoil and 44 cells normal to the airfoil (216x44 grid), selected after obtaining good agreement in the surface pressure distribution and aerodynamic coefficients at transonic conditions (Mach number of 0.729, Reynolds number of 6.5x10 6 and an angle of attack of 2.31°) as reported in Cook et al. (1979) and Slater (2002).
Optimization algorithm: A Genetic Algorithm (GA) is used to perform the optimization in this study. GA searches from multiple points in the design space, instead of moving from a single point like gradientbased methods do making it less prone to being trapped by local optima. This makes it particularly suitable for aerodynamic optimization where the function landscape is often multimodal and nonlinear because the flow field is governed by a system of nonlinear partial differential equations (Oyama et al., 2000). Furthermore, GA works on function evaluations alone and does not require derivatives or gradients of the objective function. These features make it a robust global optimization algorithm.
The GA implementation in MATLAB and its optimization toolbox is used to perform the optimization.

Sample plan:
As with all surrogate-based methods, to approximate a function f we start with a set of sample data-computed at a set of points in the domain of interest determined by a sampling plan. Selection of the sample points is a very important step towards creating a good surrogate model. If the sample plan model is too sparse or does not contain the interesting features of the design space, the resulting model will fail to resolve desirable global features. In order to improve the surrogate it is necessary to incrementally add points in an intelligent way such that the generated surface converges toward the true surface.
When dealing with large, complex design spaces, it is often unclear how many points may be necessary to resolve key features with a response surface. To get the maximum amount of information out of a minimum number of points with no a priori knowledge of the design space requires uniform sampling.
Several sampling methods are capable of producing relatively uniform samples, e.g. Latin Hypercube sampling (Forrester et al., 2008), however not many methods allow incremental uniform sampling. Latin Hypercube sampling requires a priori knowledge of how many points are desired in order to divide the domain into the appropriate number of hypercubes. This method creates nearly uniform point distributions, but requires a completely new set of data if additional points are desired. Another sampling method known as the Sobol Sequence (Sobol, 1967) has good space filling properties and allows incremental uniform sampling (Nelson et al., 2007;Forrester et al., 2008). This method is, thus, adopted in this study to create the sample plan. Since the number of sample points required to obtain an accurate surrogate is generally unknown a priori, an initial sample of 10n var (where n var is the number of design variables) is created following the suggestion of Jones et al. (1998).

SURROGATE MODEL
The surrogate model used in the study is Krigingan approximation technique that has received much attention in recent years (Sacks et al., 1989;Jones et al., 1998;Keane, 2003;Leary et al., 2003;Gano et al., 2005;Forrester et al., 2006;Huang et al., 2006;Forrester et al., 2007;Nelson et al., 2007) named after the pioneering work of the South African mining engineer D.G. Krige and introduced in engineering design work after the seminal paper by Sacks et al. (1989).
For brevity, the Kriging equations are not mentioned here. Forrester et al. (2008) contains more detail about Kriging.
The Kriging model in this study uses a constant regression term and a Gaussian correlation model. The p term is 2 for all dimensions and θ is optimized in the range 10 2 ≤θ i ≤200, i = 1, 2,…,n var . The surrogate model is created using the surrogates toolbox (Viana, 2010) for MATLAB.

Design variables:
The design process begins with an initial airfoil. The airfoil geometry is then modified by adding smooth perturbations in the form of the Hicks-Henne bump functions (Hicks and Henne, 1978). Defining normalized airfoil coordinates ψ = x/c, ζ = z/c and c as the chord length, the Hicks-Henne shape function modifies airfoil geometry by adding a linear combination of shape functions, F i and weighting coefficients, x i as follows: (1) Here t 1 locates the maximum point of the bump and t 2 controls the width of the bump. The design variables are the coefficients x i multiplying the various Hicks-Henne bump functions.
In this study, 7 bump functions are used for the upper and lower surface of the airfoil, resulting in a total of 14 design variables. The points t 1 are linearly spaced between 0 and 0.94. The range is terminated a little before the trailing edge to prevent the upper and lower edges from crossing each other and creating unrealistic geometries. A value of 10 is used for t 2 following Castonguay and Nadarajah (2007) recommendation.
To prevent large changes to the geometry, upper and lower bounds are set on x i . These are: (2) Fitness function and constraints: This is a single objective optimization problem. The airfoil is optimized for a transonic Mach number of 0.729 and a fully turbulent Reynolds number of 6.5x10 6 . The objective is to minimize the drag, c d at a constant lift, c l =0.6±0.01, subject to thickness and pitching moment constraints. RAE2822 is used as the initial airfoil to start the optimization. The optimization problem is stated as follows: Minimize: Constant lift is maintained by varying the angle of attack, α, of the airfoil. At low to moderate values of α, i.e., before flow separation and stall, c l varies linearly with α. During each evaluation of the fitness function, the airfoil is first simulated at an initially guessed α 1 and at α 2 . A third simulation at α 3 , found from a linear interpolation through (α 1 , c l1 ) and (α 2 , c l2 ), is then used to attain the desired lift. This procedure effectively removes the c l constraint from the optimization by making it a condition which each simulation must satisfy. The optimization is simplified both by the removal of the constraint and by reducing the dimensionality of the problem, since α is not used as an optimization variable. Other constraints are imposed by adding penalty terms to the fitness function: where (t/c) * is the minimum allowable thickness of 12.11% (the (t/c) max of the initial airfoil) and is -0.122 The VFO framework: A variable fidelity prediction predicts the function response in 2 steps: • The low fidelity function is evaluated to obtain an estimate • The estimate is corrected by performing a high fidelity simulation to obtain a better estimate to the function and apply a correction to the low fidelity prediction. The correction between the high and low fidelity functions is modeled by a surrogate model developed by sampling the function at a few points. To evaluate the goodness of fit over the entire function landscape, the RMSE and R 2 of the surrogate is calculated for a validation dataset generated using a uniform Latin Hypercube sample consisting of n test = n var *10 points.  Figure 1 shows the flowchart of the VFO framework used in this study. If the initial number of points used to construct the surrogate do not result in an accurate surrogate, the surrogate is updated by adding more points until some accuracy criteria is satisfied. The surrogate is then directly coupled to the GA to determine the optimum.
It has been reported that when a surrogate is used for fitness evaluation, it is very likely that the evolutionary algorithm will converge to a false optimum (Jin, 2005). A false optimum is an optimum of the approximate model, which is not one of the original fitness function. To avoid this, the VF optimum point is evaluated by the HF solver every 10 generations. If the VF lift coefficient, c l , agrees with the HF c l prediction within a 0.01 tolerance and the drag coefficient within 0.0001, the optimization is continued. Otherwise, the surrogate model is updated by performing HF evaluations on the GA population and refitting the surrogate before continuing the optimization. high and low fidelity solvers. Experimental data for the selected case is also shown for reference (Zahir and Gao, 2011) Optimization is terminated once the cumulative change in the fitness function value over 50 generations is less than or equal to 1.0 6 , or when 100 GA generations are completed (whichever comes sooner).
Parallelization: MATLAB's distributed computing toolbox provides a scheduler to distribute multiple tasks on several cores of the same computer, or on different computers. Although the CFD code used is a serial code, the unique individuals of each population are evaluated in parallel on several computers (and cores) achieving significant time savings. A heterogeneous cluster of 7 dual-core IBM-PC computers clocked at 2.5 GHz and 3.0 GHz (resulting in 14 nodes) and a head node are used in this study.

RESULTS AND DISCUSSION
Use of Euler and RANS solvers with several grid sizes was investigated previously by Zahir and Gao (2011) for use as the low-fidelity solvers. The surface pressure distribution on the RAE-2822 at Mach number of 0.729, Reynolds number of 6.5x10 6 and an angle of attack of 2.31° are shown in Fig. 2 (Zahir and Gao, 2011). Results of the HF solver (NS 216x44) are also shown. All Navier Stokes solvers use the K-ω turbulence model and an initial grid line spacing of 2x10 6 units from the airfoil surface. The Euler solvers use an initial grid line spacing of 0.01 units from the airfoil surface. Difference between the aerodynamic coefficients of the low-fidelity solvers and the HF solver is shown in Fig. 3. The computation time for one run is also shown (Zahir and Gao, 2011). Three LF solvers are selected for the VF surrogate in this study to cover the entire spectrum from slow to fast and less accurate to more accurate: Euler 60x20, Euler 160x24, RANS 160x24. Figure 4 shows the  Table 1: Aerodynamic coefficients of the VF optimized airfoils using several LF Solvers. High fidelity evaluations at the VF optimum points are also shown along with total optimization time. Results of direct optimization using the HF solver alone are also shown for reference Aerodynamic coefficients of optimized airfoils Optimization time (hours) c l c d c m K Serial Run (Zahir and Gao, 2011) 0.600 0.00643 -0.123 93.45 Parallel Run (Zahir and Gao, 2012)  It is clear from Table 1 that the RANS 160x24 solver yields the best lift-to-drag ratio, k, within 1% of Fig. 6: Geometry of the optimum airfoils. Airfoil produced by direct optimization using the HF solver, and the baseline RAE-2822 is also shown the direct optimization result. Other solvers perform well too; yielding k within 7% of the direct optimization result, while still providing a significant saving in computation time. All solvers also satisfy the thickness constraint (t/c) max ≥ 12.11%. The benefits of parallel computation are shown in Table 2 for a similar optimization problem, but without the pitching moment constraint, previously reported in Zahir and Gao (2011) and Zahir and Gao (2012). The parallel computations are 24 times faster than the serial computations a very significant improvement.
The surface pressure distributions on the optimum airfoils are shown in Fig. 5 and the airfoil geometries are shown in Fig. 6. It is seen that the optimum pressure distributions and airfoil geometries produced using different LF solvers are quite similar. This is also reflected in the aerodynamic coefficients shown previously in Table 1.

CONCLUDING REMARKS
The VFO method used in this study is, both, efficient and capable of finding the optimum that closely agrees with the results of high-fidelity optimization alone. Modification of the authors' original VF framework to incorporate parallel evaluation yields significant time savings. Euler and Navier-Stokes solvers evaluated on low resolutions grids are good candidates for LF solvers and produce results close to the HF optimum with significant time savings.