Aerodynamic optimization of an airfoil using gradient based method

AERODYNAMIC OPTIMIZATION OF AN AIRFOIL USING

GRADIENT BASED METHOD

Masoud Mirzaei†, Jafar Roshanian*, Seyeded Nasrin Hosseini**

†

Assistant Professor, Aerospace Faculty, K.N.Toosi University of Technology, Tehran pars 4th square, Vafadar Ave., Tehran, Iran

e-mail: mirzaei@kntu.ac.ir

Web page: http://sahand.kntu.ac.ir/~mirzaei/

*

Assistant Professor, MDO Lab, K.N.Toosi University of Technology Tehran pars 4th square, Vafadar Ave., Tehran, Iran

e-mail: roshanian@kntu.ac.ir

**

Graduate Student, Aerospace Faculty, K.N.Toosi University of Technology, Tehran pars 4th square, Vafadar Ave., Tehran, Iran

e-mail: as_nasrin@yahoo.com

Key words: CFD, Optimization, Gradient, Objective function, Design variables

Abstract. A gradient based method is presented for optimization of an airfoil configuration. The flow is governed by two dimensional, compressible Euler equations. A finite volume code based on unstructured grid is developed to solve the equations. The procedure is carried out for optimizing an airfoil with initial configuration of NACA 0012. The advantage of this technique over the other gradient based methods is its speed of converging.

1 INTRODUCTION

Computational fluid dynamics (CFD) and numerical optimization techniques are being used widely in the field of aerodynamic shape optimization. This growing interest is owing to the fact that CFD is making fast progress using advanced computer technology and available efficient numerical algorithms. This helps faster computation of the flow field which is necessary for the optimization. Also, efficient numerical algorithms for optimization are available so that the combined effort makes it possible to compute optimum solution in considerable time. This helps saving the cost occurred in experimental methods [1].

Finite difference and Newton's methods are examples of gradient methods [3]. Adjoint method is also a gradient based method which have been developed and applied in recent years. Pironneau [4] was the first one used adjoint equations for design purpose in fluid dynamics. Jameson [5] applied this method in aerodynamic problems. Other researchers such as Elliot [6], Qioa [7], Shenoy [8], Xie [9], Nadarajah [10], Zingg, Gatsis [11] implemented adjoint equations in aerodynamic design problems.

In this study finite difference method is used to estimate the gradient of the objective function. Here finding an optimized airfoil profile is required. To do this, optimization is initiated from a given airfoil and the flow field is solved; then, gradient of objective function is calculated with respect to the design variables. With new design variables the flow field is solved again, and this procedure is iterated to reach the minimum value of the objective function. The final design variables which lead to the minimum value objective function are optimum ones. Design variables are coefficients of a polynomial equation defining upper and lower airfoil profiles. Aerodynamic characteristics such as

L

C

1 _{and are assumed as an} objective function. Optimization of a given airfoil (NACA 0012) is studied using these two objective functions.

D

C

2 OPTIMIZATION APPROACH

As stated previously, a gradient based optimization method is adopted. The method is a successive quadratic programming which is generally used to solve non linear problems. Normally the optimization problem is stated as follows:

u l R x x x x x f n ≤ ≤ ∈ ) (

min

f is assumed to be continuously differentiable. The method is based on iterative formulation and solution of quadratic programming of the related problems. It obtains the related problems using a quadratic approximation of the Lagrangian. That is:

k u k l T k k T R x x x d x x d x f d B d n − ≤ ≤ − ∇ + ∈ ) ( 2 1

min

Where is a positive definite approximation of the Hessian and is the current iteration. Let be the solution of the problem. A line search is used to find a new point :

k B xk k d xk+1 ] 1 , 0 ( , 1 = + ∈ + k λ k λ k x d x (3)

point. Here, the augmented Lagrange function is used as the objective function. When optimality is not achieved, B_k is updated.

This approach needs initial data, design variables and objective function. NACA0012 profile is adopted as initial data, coefficients of the polynomial defining upper and lower profiles of the airfoil, are considered as design variables and

L C 1 _{or are assumed as} objective function. D C

3 FLOW FIELD GOVERNING EQUATIONS

Euler equations are adopted as the flow field governing equations and have essential role in determination of the objective function. These equations are developed from conservation of mass, momentum and energy laws. Vector form of these equations is:

0 = ∂ ∂ + ∂ ∂ + ∂ ∂ y F x E t U (4) Where: ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = t e v u U ρ ρ ρ ρ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + = t vh P v vu v E ρ ρ ρ ρ 2

E and Fare flux vectors in x and directions respectively. U is vector of dependent variables, and and are total energy and total enthalpy respectively.

y

t

e h_t

The above set of equations contains four equations with five unknowns, to close this set; we need an extra equation which will be the equation of state. For an idea gas this equation is:

e p ) 1 ( − = γ ρ (6)

Where is internal energy and e γ is specific heat ratio.

4 NUMERIC

Finite volume technique is adopted for spatial disertization of the governing equations. After integration from equation 4 over triangular control volumes shown in figure 1, we have:

the length of the control volume edge, and denotes area of the control volume j. The next

step is determination of fluxes on the edge of the control volume. Here the inviscid fluxes are approximated using Roe’s scheme. This scheme belongs to upwind schemes category so called flux difference splitting that uses approximate Riemann problems for flux calculations [12]. One of the distinct characteristic of this scheme is its excellent capability in capturing of flow field discontinuities such as shock waves.

j

A

After special integration from Euler equations, the equations are converted to ordinary differential equations:

(

)

Where, R is numerical flux function. Equation 8 is an ordinary differential equation which may be solved using classical methods such as 4th order Runge-Kutta method [13]:

4 ,..., 1 , 1 4 1 1 = + − Δ + = − − k R k t U Uk k k (9)

5 GRID GENERATIONS

As stated before, the discretization of the equations is carried out in triangular cells. Here Delauney triangulation method [14] is chosen to discretize the physical domain. This approach is designed for applying an optimum method in node connections in triangular unstructured mesh. The most important specification of this method is its efficiency and mesh generated quality. Since the geometry of the airfoil varies during optimization process, new cell arrangement is needed for new geometry. This changing is applied by method named spring-based smoothing method [15]. In the spring-based smoothing method, the edges between any two mesh nodes are idealized as a network of interconnected springs. The initial spacing of the edges before any boundary motion constitutes the equilibrium state of the mesh. A displacement at a given boundary node will generate a force proportional to the displacement along all the springs connected to the node. Figure 2 represent the performance of the method.

6 RESULTS

Since the flow field solution performs major role in evaluation of the objective function; the flow field solver should be validated. The results of the present solver are compared with those of [16]. Figure 3 gives the pressure coefficient distribution over lower and upper surface of NACA 0012 airfoil for free stream conditions ofM =0.7, it can be seen that the present results have good agreement with the experimental data of [16]. Figure 4 shows the pressure contours of the present solver for this case.

0 49 . 1 = α

function of L

C

1 _{. Table 1 shows the results of this optimization. In this table the initial and} optimized values of design variables and objective function are presented. It should be noted that the airfoil profile is:

) ( 2 . 0 4 5 3 4 2 3 2 1 x a x a x a x a x a t y = + + + + (10) 5 1,...,a

a are design variables, and t is for the airfoil. Figure 5 .a depicts the convergence history of the objective function. It can be seen that the method is converged after 13 times function evaluation. Comparing with the other methods it converged faster; hence, it is very time saving approach.

12 . 0 1 a a ₂ a 3 a 4 a 5 L C 1 L C Initial (NACA 0012) 0.2969 -0.126 -0.3516 0.2843 -0.1015 3.597 0.278 Optimized 0.3711 -0.1575 -0.4395 0.3553 -0.1268 3.194 0.313

Table 1: Design results (case 1)

Pressure distribution and airfoil geometry of initial and optimized profiles are shown in figure 5.b and 5.c respectively. According to the figure 5.c the optimized airfoil is thicker and consequently produces more lift.

In another test case, NACA 0012 is considered as an initial airfoil again and free stream conditions of and are chosen, the objective function is considered . Table 2 shows the results of this optimization, the initial and optimized values of design variables and objective function are presented. Figure 6.a displays the convergence history of objective function. The pressure distribution and airfoil geometry of initial and optimized profiles are shown in figure 6.b and 6.c respectively. According to figure 6.c, the optimized airfoil is thinner than the initial one; consequently, total drag is smaller and shock location has moved nearer to leading edge; in addition, its strength has been reduced too. Optimized solution has been obtained after 15 design cycles, with respect to 6 design variables and as stated before it is really a fast method than other similar methods.

8 . 0 = M 0 2 . 2 = α C_D 1 a a ₂ a 3 a 4 a 5 C D α Initial (NACA 0012) 0.2969 -0.126 -0.3516 0.2843 -0.1015 0.045 2.2 Optimized 0.2226 -.0945 -0.2637 0.2132 -0.0761 0.016 1.4

Table 2: design result (case 2)

7 CONCLUSIONS

of L

C

1 _{, optimization leads to 11% increase in lift coefficient after 13 design cycles. For the} second case with optimization function of the optimized case shows 64 % reduction in drag coefficient in comparison with initial airfoil after 15 design cycles.

D

C

8 REFERENCES

[1] Subhendu Bikash Hazra, An Efficient Method for Aerodynamic Shape Optimization, AIAA 2004-4628, 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization

Conference 30 August- 1 September 2004, Albany, New York.

[2] Gen, M., Cheng, R., “Genetic Algorithm & Engineering Optimization”, John Wily & Sons, 2000.

[3] http://en.wikipedia.org/wiki/multidisciplinary_design_optimization

[4] Pironneau, O., An Optimum Design in Fluid mechanics, J. Fluid Mech, 1974.

[5] A. Jameson, Aerodynamic design via control theory, Journal of Scientific Computing 1998.

[6] J.K. Elliot, Aerodynamic Optimization Based on the Euler and Navier-Stokes Equations Using Unstructured Grids, Ph.D. Thesis, Massachusetts Institutes of Technology, 1998 . [7] Z.D. Qiao, X.D. Yang, X.L Qin, B. Zhu, Numerical Optimization Design by Solving

Adjoint Equations, ICAS 2002 Congress.

[8] A.R. Shenoy. Optimization Techniques Exploiting Problem Structure: Application to Aerodynamic Design, Ph.D. Thesis, Virginia Polytechnic Institute and State University, 1997.

[9] L. Xie, Gradient-Based Optimum Aerodynamic Design Using Adjoint Methods, Ph.D. Thesis, Virginia Polytechnic Institute and State University, 2002.

[10] S.K. Nadarjah, the Discrete Adjoint Approach to Aerodynamic Shape Optimization, Ph.D. Thesis, Stanford University, 2003.

[11] J.Gatsis, D.W. Zingg, a Fully Coupled Algorithm for Aerodynamic Design Optimization, University of Toronto Institute for Aerospace Studies.

[12] Roe, P.L. The Use of the Riemann Problem in finite difference schemes, in Proc. 8th International Conference on Numerical Methods in Fluid Dynamics, Berlin, Lecture Notes in Physics, vol. 141, Springer-Verlag, New York , pp.354-9, 1981.

[13] Mirzaei M., Shahverdi M. A comparison between Kinetic Flux Vector Splitting and Flux Difference Splitting Methods in Solution of Euler Equations, CFD2004 conference, Canada, Ottawa, May 2004

[14] Klaus A.Hoffmann, Steve T.Chiang, Computational Fluid Dynamics for Engineers, volume II.

[15] Farhat, C., Degand, C., Koobus, B., and Lesoinne, M., “An Improved method of Spring Analogy for Dynamic Unstructured Fluid Meshes”, AIAA Paper 98-2070, 1998.

[16] http://xoptimum.narod.ru/results/compressible/naca0012

[17] Marian Nemec, David Zingg, “Multi-Point and Multi-Objective Aerodynamic Shape Optimization”, 9th AIAA, Sept. 4-6, 2002/Atlanta, Georgia.

Lectures at Von Karman Institute, Brussels February 6, 2003, Stanford, California.

[19] Eric J.Nielson, W.Kyle Anderson, Recent Improvement in Aerodynamic Design Optimization on Unstructured Meshes, AIAA Journal, Vol.40, No.6, 2002, pp.1155-1163.

[20] Chasidic Leoviriyakit, Antony Jameson, “Aerodynamic Shape Optimization of Wings Including Planform variations”, 41st Aerospace Sciences Meeting and exhibit, Jan 6-9 2003/Reno, Nevada.

[21] M.Nemec, D.W.Zingg, “Newton-Krylov Algorithm for Aerodynamic Design Using the Navier –Stokes Equations”, AIAA Journal, Vol.40, No.6, 2002, pp.1146-1154

[22] A.Jameson, S.Kim, S.Shankaran, K. Leoviriyakit, “Aerodynamic Shape Optimization: Exploring the Limits of Design”, KSAS 1st International Sessions in 2003 Fall Conference, Nov 14-15, 2003, GyeongJu, Korea.

X Y 0 0.25 0.5 0.75 1 -0.4 -0.2 0 0.2 0.4

Figure 1. Triangular control volume

X Y 0 0.025 0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04

x -C p 0 0.25 0.5 0.75 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 1.25 _Experiment Solver

Figure 3: Comparison of experimental Cp and Calculated Cp by solver X Y 0 0.5 1 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Figure 4: Pressure contours of the present solver for validation case Design Cycle O b je ct iv e F u n ct io n (1 /C L ) 5 10 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.55 3.6

Figure 5.a: History of Objective Function (case 1)

Design cycle O b je ct ive F un ct io n (CD ) 5 10 15 0.015 0.02 0.025 0.03 0.035 0.04 0.045

x -Cp 0 0.25 0.5 0.75 1 -1 -0.5 0 0.5 1 Optimized NACA 0012

Figure 5.b: Pressure Distribution (case 1)

x -C p 0 0.25 0.5 0.75 1 -1 -0.5 0 0.5 1 Optimized NACA 0012

Figure 6.b: Pressure Distribution (case 2)

x y -0.25 0 0.25 0.5 0.75 1 1.25 -0.3 -0.2 -0.1 0 0.1 0.2 Optimized NACA 0012

Figure 5.c: Airfoil Configurations (case 1)