Optimization of Airfoils for Minimum Pitching Moment and Compressibility Drag Coefficients

ABSTRACT:

This paper concerns a numerical optimization method for designing airfoils based on adjoint method. The goal of present work is to reduce the compressibility drag or pitching moment of transonic airfoils without compromising on the lift coefficient. A new cost function based on this requirement is defined and the corresponding adjoint equations are discussed in details. At the end, by demonstrating some numerical results, we show that this technique is capable of converging to the optimum design point corresponding to the initial geometry of the airfoil.

KEYWORDS:
Optimization; Pitching moment; Transonic airfoils; Adjoint equations; Euler equations

INTRODUCTION

Improvement of aerodynamic performance of helicopter blade designs has been one of the most important research areas in rotorcraft aerodynamics. Low-pitching-moment airfoils found application primarily as helicopter rotor blades, but some attention has been given to the advantages of low pitching-moment sections for a "span-loader" vehicle. For such applications, a symmetric airfoil could conceivably be employed, but cambered airfoils can offer significant advantages (Barger, 1975). Therefore, designing cambered airfoils with low pitching-moment is attractive in helicopter blade design to achieve small control forces in rotor controls, and this task can be conducted via optimizing procedures.

Also, shock waves in transonic flow not only increase the wave drag but also cause unfavourable flutter or buffet to fixed wings. Helicopter rotor blades in forward flight are also frequently exposed to strong unsteady shock waves at the blade tip region. These shock waves increase the required torque, and become a source of undesirable noise and vibration. Thus, elimination or possible reduction of the strength of these shock waves would be desirable for the enhancement of the performance of helicopters and airplanes.

In order to improve airfoil performance under different flight conditions and to make the performance insensitive to off-design condition at the same time, Yu et al. (2010) have developed a multi-objective optimization approach considering robust design and applied it to airfoil tailoring. They used non-uniform rational B-spline (NURBS) representation, control points and related weights around airfoil as design variables. They could obtain a set of non-dominated airfoil solutions with robustness, by adopting multi-objective genetic algorithm that is based on non-dominated sorting. They showed that the multi-objective robust design method makes the airfoil performance robust under different off-design conditions.

In the gradient-based optimization design, the gradients of cost/constraint functionals with respect to design variables are important information in the design processes (Xie, 2002). A traditional approach to calculate gradients is the finite difference method which involves finite differencing performance functionals computed from the solutions to the governing equations with perturbed parameter values. The Finite difference method (FDM) with complex variables was explored by Anderson et al. (2001) to improve its accuracy. Even though it is very easy to implement FDM in program coding, its prohibitive computational costs (many times solving governing equations) motivate other time-efficient methods to calculate reduced gradients.

A way of performing an efficient aerodynamic shape optimization is to regard the design problem as a control problem in which the control is the shape of the boundary. This approach to optimal aerodynamic design was introduced by Jameson (1988, 1990) who examined the design problem for compressible flow with shock waves, and devised adjoint equations to determine the gradient for both potential flow and also flows governed by the Euler equations (Jameson, 1995; Reuther and Jameson, 1994).

Also, particular interest has been given to adjoint methods (Reuther et al., 2001) in which the gradient, regarding an arbitrarily large number of parameters, can be calculated with roughly the same computational cost as two flow solutions. Once the gradient has been calculated, a descent method can be used to determine a shape change which will make an improvement in the design.

Based on the idea of adjoint method, an optimum aerodynamic design technique is presented by Ying et al. (2011), which can be applied to the optimum problems with a large number of design variables. The key of this method lies in that the optimization process is regarded as an unsteady evolution, i.e., the optimization is executed, simultaneously with solving the unsteady flow governing equations and adjoint equations. During the unsteady evolution, the airfoil surface is moving with time, and the dynamic grid technique is introduced in the grid generation to improve the computational efficiency (Ying et al., 2011).

Nadarajah and Tatossian (2008) also have presented an adjoint method for the multi-objective aerodynamic shape optimization of unsteady viscous flows. Their method is beneficial when both the pitching angle and free stream Mach number are sinusoidally varied. They showed that the multi-objective cost function was able to preserve both the lift and pitching moment while simultaneously decreasing the overall drag (Nadarajah and Tatossian, 2008).

Gradient-based techniques are efficient direct aerodynamic shape optimization (ASO) methods for both incompressible and compressible flows, and many of them use continuous adjoint methods.

More recently, the introduction of surrogate-based optimization (SBO) methods to ASO have been successful in reducing the total computational cost. The overall objective of using SBO methods is to reduce the number of evaluations of the high-fidelity models, and thereby making the optimization process more efficient. Leifsson and Koziel (2010) introduced a computational design methodology which exploits surrogates constructed using low-fidelity flow analysis models and shape-preserving response prediction technique and demonstrated that their approach allows a rapid design improvement of airfoils at a very low computational cost corresponding to a few evaluations of the high-fidelity model (Leifsson and Koziel, 2010).

Then they replaced the direct optimization of an accurate, high fidelity airfoil model by an iterative re-optimization of a corrected low-fidelity model (Leifsson et al., 2011). Their low-fidelity model is based on the same governing fluid flow equations as the high-fidelity one, but uses coarser discretization and related convergence criteria. With these techniques, they succeeded to improve the overall robustness of the optimization algorithm.

ASO problems in compressible viscous flow were also performed using simultaneous pseudo-time stepping (Hazra et al., 2005). Hazra et al. (2005) could optimize airfoil surface by use of a preconditioner for convergence acceleration which stems from the reduced sequential quadratic programming (SQP) methods.

This paper describes the implementation of adjoint technique for designing transonic airfoils with minimum compressible drag or pitching moment coefficients while preserving the initial lift coefficient. For this purpose, the Euler and adjoint equations are solved sequentially on a two-block structured grid about the section. The finite-volume scheme of Jameson et al.(1981) is used to discrete governing equations. This scheme employs a cell centered discretization technique. A multistage time-stepping algorithm is used to advance the solution in time. A key feature in this scheme is applying the dissipative terms at the end of each time step for solving the flow and adjoint equations. The magnitude of the dissipative terms is adapted to the local properties of the flow by means of a Jameson sensor (Nadarajah and Jameson, 2001), based on the local pressure gradient. Acceleration techniques are applied to obtain faster steady-state convergence. These methods include local time stepping, variable coefficient implicit residual smoothing, and multi grid technique.

GOVERNING EQUATIONS

Let (u,v) denote the velocity components in the Cartesian coordinate system (x,y). Compressible Euler equations can be formulated as:

(1) $\frac{\mathit{\partial }\mathit{W}}{\mathit{\partial }\mathit{l}}+\frac{\mathit{\partial }{\mathit{f}}_{\mathit{X}}}{\mathit{\partial }\mathit{x}}+\frac{\mathit{\partial }{\mathit{f}}_{\mathit{y}}}{\mathit{\partial }\mathit{y}}=0$

Where:

(2) $\mathit{W}=\left\{\begin{array}{c}\mathit{\rho }\\ \mathit{\rho }\mathit{u}\\ \mathit{\rho }\mathit{v}\\ \mathit{\rho }\mathit{e}\end{array}\right\},{\mathit{f}}_{\mathit{X}}=\left\{\begin{array}{c}\mathit{\rho }\mathit{u}\\ \mathit{\rho }{\mathit{u}}^2+\mathit{p}\\ \mathit{\rho }\mathit{uv}\\ \mathit{\rho }\mathit{uh}\end{array}\right\},{\mathit{f}}_{1\prime }=\left\{\begin{array}{c}\mathit{\rho }\mathit{v}\\ \mathit{\rho }\mathit{vu}\\ \mathit{\rho }{\mathit{v}}^2+\mathit{p}\\ \mathit{\rho }\mathit{vh}\end{array}\right\}$

Defining the flux vecto r as ƒ=ƒ_xî+ƒ_yĵ, where (î,ĵ) are the unit vectors in the (x, y) coordinate system, Eq. 1 can be written as:

(3) $\frac{\mathit{\partial }\mathit{W}}{\mathit{\partial }\mathit{t}}+\nabla .\mathit{f}=0$

Equation 3 can be solved by the finite-volume method described in the proceeding section.

NUMERICAL METHOD

In order to apply the finite-volume method, Eq. 3 is written in the integral form as:

(4) $\frac{\mathit{\partial }}{\mathit{\partial }\mathit{t}}\underset{\mathit{V}}{\int }\mathbf{W}\mathit{dV}+\underset{\mathit{\partial }}{\int }\underset{\mathit{S}}{\int }\hat{\mathit{n}}.\mathit{f}\mathit{ds}+\underset{\mathit{\partial }}{\int }\underset{\mathit{S}}{\int }{\mathit{u}}_{\Omega }.\hat{\mathit{n}}\mathit{W}\mathit{ds}=\int \mathit{S}\mathit{dV}$

The computational domain is meshed by a multi-block structured quadrilateral grid. Applying Eq. 4 to a typical (i,j) cell, for a cell center algorithm we obtain a system of ordinary differential equations as:

(5) $\frac{\mathit{d}}{\mathit{dt}}({\forall }_{\mathit{ij}}{\mathit{W}}_{\mathit{ij}})+F_{\mathit{ij}}=0$

Where ∀_ij is the cell area and F_ij represents the net convectional fluxes out of each cell.

The convective fluxes can be determined by a second order central scheme. However, the resultant scheme is not dissipative, and therefore, undamped oscillations at odd and even mesh points can be developed during the computation. To suppress the tendency for odd-even decoupling and to prevent the appearance of oscillations in regions containing severe pressure gradients near shock waves and stagnation points, the finite-volume scheme is modified by the addition of artificial dissipative terms as below:

(6) $\frac{\mathit{d}}{\mathit{dt}}({\forall }_{\mathit{ij}}{\mathit{W}}_{\mathit{ij}})+F_{\mathit{ij}}-{\mathit{D}}_{\mathit{ij}}=0$

Where D_ij denotes the dissipative fluxes. Jameson has established that an effective form of dissipative terms for flows with discontinuities is a blend of second and fourth differences with coefficients which depend on the local pressure gradient (Jameson et al., 1981). Equation 6 can be integrated in time by a multistage Runge-Kutta scheme.

BOUNDARY CONDITIONS

Special care must be given to the boundary conditions when solving Euler equations on a multi-block grid. When using more than a single block for a two-dimensional body, the attachment boundary condition should be used for the neighboring cells at the block interfaces. This can be done by considering a layer of ghost cells around each block. The boundary condition at the airfoil surface is the tangency condition of flow, that is:

(7) $\overrightarrow{V}.\hat{n}=0$

Where V⃗ is the velocity vector and n̂ is the unit vector normal to the airfoil surface.

Far-field boundary condition is also applied at the outer surface of the domain. The treatment of the far-field boundary condition is based on the concept of Riemann invariants for a one-dimensional flow normal to the boundary (Jameson and Baker, 1983).

DERIVATION OF THE INVISCID ADJOINT TERMS

The main objective of the present study is to achieve a modified airfoil which minimizes the inviscid compressibility drag without compromising the pitching moment under the constraint of maintaining the desired lift level. For this purpose, using the Lagrange multiplier, the cost function is defined as:

(8) $\mathit{I}={\mathit{k}}_1\mathit{Cd}+{\mathit{k}}_2{\mathrm{Cm}}^2+\sigma (\mathit{Cl}-{\mathit{Cl}}_{\mathit{d}})$

Here k₁ and k₂ are the weight coefficients by which the airfoil shape can be modified for a particular application. When k₁ is chosen to be much larger than k₂, a reduction in the drag coefficient is the main priority of the designer while using a larger k₂ leads to an airfoil with smaller pitching moment coefficient. In this equation Cl_d indicates the desired value of Cl. The variation of cost function in relation to the angle of attack and other parameters can be written as:

(9) $\begin{array}{c}\frac{\partial \mathit{I}}{\partial \mathit{\alpha }}={\mathit{k}}_1\frac{\partial \mathit{Cd}}{\partial \mathit{\alpha }}{\mid }_{\mathit{G}}+{\mathit{k}}_2\mathit{Cm}\frac{\partial \mathit{Cm}}{\partial \mathit{\alpha }}{\mid }_{\mathit{G}}+\mathit{\sigma }\frac{\partial \mathit{Cl}}{\partial \mathit{\alpha }}{\mid }_{\mathit{G}}=0\hfill \\ \frac{\partial \mathit{I}}{\partial \mathbf{G}}{\mid }_{\mathit{Cl}}={\mathit{k}}_1\frac{\partial \mathit{Cd}}{\partial \mathbf{G}}{\mid }_{\mathit{\alpha }}+{\mathit{k}}_2\mathit{Cm}\frac{\partial \mathit{Cm}}{\partial \mathbf{G}}{\mid }_{\mathit{\alpha }}+\mathit{\sigma }\frac{\partial \mathit{Cl}}{\partial \mathbf{G}}{\mid }_{\mathit{\alpha }}\hfill \end{array}$

Where |_(.) denotes the partial derivative at the frozen quantity (.) and G denotes the flow or design variables. From Eq. 9, the Lagrange multiplier σ can be obtained as:

(10) $\sigma =-{\mathit{k}}_1\frac{\mathit{\partial }\mathit{Cd}\mathit{/}\mathit{\partial }\mathit{\alpha }{|}_G}{\mathit{\partial }\mathit{Cl}\mathit{/}\mathit{\partial }\mathit{\alpha }{|}_G}-{\mathit{k}}_2\frac{\mathit{\partial }\mathit{Cm}\mathit{/}\mathit{\partial }\mathit{\alpha }{|}_G}{\mathit{\partial }\mathit{Cl}\mathit{/}\mathit{\partial }\mathit{\alpha }{|}_G}$

Substituting into Eq. 9 yields:

(11) $\frac{\mathit{\partial }\mathit{I}}{\mathit{\partial }\mathbf{G}}{|}_{\mathit{Cl}}={\mathit{k}}_1\left(\frac{\mathit{\partial }\mathit{Cd}}{\mathit{\partial }\mathbf{G}}{|}_{\mathit{\alpha }}+{\Omega }_1\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathbf{G}}{|}_{\mathit{\alpha }}\right)+{\mathit{k}}_{2}C\mathit{m}\left(\frac{\mathit{\partial }\mathit{Cm}}{\mathit{\partial }\mathbf{G}}{|}_{\mathit{a}}+{\Omega }_2\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathbf{G}}{|}_{\mathit{\alpha }}\right)$

Where: ${\Omega }_1=-\frac{\mathit{\partial }\mathit{Cd}\mathit{/}\mathit{\partial }\mathit{\alpha }{\mathit{|}}_G}{\mathit{\partial }\mathit{Cl}\mathit{/}\mathit{\partial }\mathit{\alpha }{\mathit{|}}_G}\mathit{,}{\Omega }_2=-\frac{\mathit{\partial }\mathit{Cm}/\mathit{\partial }\mathit{\alpha }{|}_G}{\mathit{\partial }\mathit{Cl}/\mathit{\partial }\mathit{\alpha }{|}_G}$

The airfoil lift, drag and pitching moment coefficients can also be written as:

(12) $\begin{array}{c}\mathit{Cl}=-\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}{\mathit{S}}_{22}\mathit{Cos}\mathit{\alpha }+\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}{\mathit{S}}_{21}\mathit{Sin}\mathit{\alpha }\hfill \\ \mathit{Cd}=-\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}{\mathit{S}}_{21}\mathit{Cos}\mathit{\alpha }-\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{P}}{\mathit{S}}_{22}\mathit{Sin}\mathit{\alpha }\hfill \\ \mathit{Cm}=-\frac{1}{\mathit{c}}\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}(\mathit{x}-{\mathit{x}}_{\mathit{ref}}){\mathit{S}}_{22}-\frac{1}{\mathit{c}}\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}{\mathit{yS}}_{21}\hfill \end{array}$

Where x_ref is the reference point for computing the pitching moment and S₂₁ and S₂₂ are the components of unit vector normal to the airfoil surface.

Taking variation from the steady state flow equations in the computational domain leads to:

(13) $\frac{\partial }{\partial {\mathit{\xi }}_{\mathit{i}}}\mathit{\delta }{\mathit{F}}_{\mathit{i}}=0$

Where: $\mathit{\delta }{\mathit{F}}_{\mathit{i}}={\mathit{C}}_{\mathit{i}}\mathit{\delta }\mathit{w}+\mathit{\delta }{\mathit{S}}_{\mathit{ij}}{\mathit{f}}_{\mathit{j}},{\mathit{C}}_{\mathit{i}}={\mathit{S}}_{\mathit{ij}}\mathit{,}{\mathit{A}}_{\mathit{j}},{\mathit{A}}_{\mathit{j}}=\frac{\mathit{\partial }{\mathit{f}}_{\mathit{j}}}{\mathit{\partial }\mathit{w}}$

Multiplying Eq. 13 by a vector of co-state variables ψ and integrating over the domain,

(14) $\underset{D}{\int }{\mathit{\psi }}^{\mathit{T}}\left(\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{\xi }}_{\mathit{k}}}\mathit{\delta }{\mathit{F}}_{\mathit{k}}\right)\mathit{d}{D}_{\xi }=0$

Finally, this equation can be integrated by parts to give:

(15) $\underset{D}{\int }\left(\frac{\mathit{\partial }{\mathit{\psi }}^{\mathit{T}}}{\mathit{\partial }{\mathit{\xi }}_{\mathit{k}}}\mathit{\delta }{\mathit{F}}_{\mathit{k}}\right)\mathit{d}{\mathbf{D}}_{\xi }=\underset{\mathit{B}}{\int }{\mathit{\psi }}^{\mathit{T}}(\mathit{\delta }{\mathit{F}}_{\mathit{k}}.{\mathit{n}}_{\mathit{k}})\mathit{d}\mathbf{E}$

Thus, the variation in the cost function can be written as:

(16) $\begin{array}{c}\mathit{\delta }\mathit{I}={\mathit{k}}_1\{-\underset{\mathit{B}}{\int }\delta {\mathit{C}}_{\mathit{p}}{\mathit{S}}_{21}\mathit{Cos}\mathit{\alpha }-\underset{\mathit{B}}{\int }{\mathit{C}}_{\mathit{p}}\delta {\mathit{S}}_{21}\mathit{Cos}\mathit{\alpha }-\hfill \\ \int \delta {\mathit{C}}_{\mathit{p}}{\mathit{S}}_{22}\mathit{Sin}\mathit{\alpha }-\int {\mathit{C}}_{\mathit{p}}\delta {\mathit{S}}_{22}\mathit{Sin}\mathit{\alpha }{\}}_{+}{\mathit{k}}_1{\Omega }_1\hfill \\ \{-\int \delta {\mathit{C}}_{\mathit{p}}{\mathit{S}}_{22}\mathit{Cos}\mathit{\alpha }-\underset{B}{\int }{\mathit{C}}_{\mathit{p}}\delta {\mathit{S}}_{22}\mathit{Cos}\mathit{\alpha }+\hfill \\ \underset{B}{\int }\mathit{\delta }{\mathit{C}}_{\mathit{p}}{\mathit{S}}_{21}\mathit{Sin}\mathit{\alpha }+\underset{B}{\int }{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{S}}_{21}\mathit{Sin}\}\mathit{\alpha }+\frac{{\mathit{k}}_2}{\mathit{c}}\mathit{Cm}\underset{\mathit{B}}{\int }[\mathit{\delta }{\mathit{C}}_{\mathit{p}}\hfill \\ (\mathit{x}-{\mathit{x}}_{\mathit{rcf}}){\mathit{S}}_{22}+{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{xS}}_{22}+{\mathit{C}}_{\mathit{p}}(\mathit{x}-{\mathit{x}}_{\mathit{rcf}})\mathit{\delta }{\mathit{S}}_{22}-\hfill \\ \mathit{\delta }{\mathit{C}}_{\mathit{p}}{\mathit{yS}}_{21}+{\mathit{k}}_2{\Omega }_2\mathit{Cm}\{\underset{\mathit{B}}{\int }[\mathit{\delta }{\mathit{C}}_{\mathit{p}}\mathit{Sin}\mathit{a}{\mathit{S}}_{21}+\hfill \\ {\mathit{C}}_{\mathit{p}}\mathit{Sin}\mathit{\alpha }\mathit{\delta }{\mathit{S}}_{21}-\mathit{\delta }{\mathit{C}}_{\mathit{p}}\mathit{Cos}\mathit{\alpha }{\mathit{S}}_{22}-{\mathit{C}}_{\mathit{p}}\mathit{Cos}\mathit{\alpha }\mathit{\delta }{\mathit{S}}_{22}\mathrm{]\}}\hfill \\ +\underset{B}{\int }{\mathit{\psi }}^{\mathit{T}}(\mathit{\delta }{\mathit{F}}_{\mathit{k}}.{\mathit{n}}_{\mathit{k}})\mathit{d}B-\underset{D}{\int }\frac{\partial {\mathit{\psi }}^{\mathit{T}}}{\partial {\mathit{\xi }}_{\mathit{k}}}\{{\mathit{C}}_{\mathit{k}}\mathit{\delta }\mathit{w}+\mathit{\delta }{\mathit{S}}_{\mathit{kl}}{\mathit{f}}_{\mathit{l}}\}\mathit{d}{\mathbf{D}}_{\mathit{\xi }}\hfill \end{array}$

and n₁=0 on boundary B.

Eliminating the term which contains δw from the above integrals, the adjoint equation can be obtained as:

(17) $\frac{\mathit{\partial }\mathit{\psi }}{\mathit{\partial }\mathit{t}}-{\mathit{C}}_{\mathit{k}}^{\mathit{T}}.\frac{\mathit{\partial }\mathit{\psi }}{\mathit{\partial }{\mathit{\xi }}_{\mathit{k}}}=0$

Because n₁=0 on airfoil surface, the only remaining component of F_k is:

(18) ${\mathit{F}}_2={\mathit{S}}_{21}{\mathit{f}}_1+{\mathit{S}}_{22}\mathit{f}$

Where from the tangent boundary condition of flow on the airfoil surface, the variation of F₂ becomes:

(19) $\mathit{\delta }{\mathit{F}}_2=\left\{\begin{array}{c}0\\ \delta S_{21}\\ \delta {\mathit{S}}_{22}\\ 0\end{array}\right\}\mathit{p}+\left\{\begin{array}{c}0\\ S_{21}\\ {\mathit{S}}_{22}\\ 0\end{array}\right\}\delta \mathit{p}$

If the coordinate transformation is such that δS is negligible in the far field, then the boundary condition of adjoint Eq. 17 can be obtained as:

(20) $\begin{array}{c}{S}_{21}{\mathit{\psi }}_2+S_{22}{\mathit{\psi }}_3+\frac{{\mathit{k}}_1}{{\mathit{q}}_{\infty }}\mathrm{\{(}{S}_{2I}\mathit{Cos}\mathit{\alpha }+S_{22}\mathit{Sin}\mathit{\alpha })\hfill \\ +{\Omega }_1(S_{22}\mathit{Cos}\mathit{\alpha }-S_{21}\mathit{Sin}\mathit{\alpha })+\frac{{\mathit{k}}_2}{{\mathit{q}}_{\infty }\mathit{c}}\mathrm{Cm}\hfill \\ \{\frac{1}{\mathit{c}}[y{S}_{21}-(\mathit{x}-{\mathit{x}}_{\mathit{rcf}})S_{22}]+{\Omega }_2\mathrm{\}(}{S}_{22}\mathit{Cos}\mathit{\alpha }-\hfill \\ S_{21}\mathit{Sin}\mathit{\alpha })=0\hfill \end{array}$

Where: ${\mathit{q}}_{\infty }=\frac{\mathit{\gamma }}{2}\mathit{M}_{\infty }{}^2{\mathit{P}}_{\infty }.$. Finally we have:

(21) $\begin{array}{c}\mathit{\delta }\mathit{I}=-{\mathit{k}}_1\{\underset{\mathbf{B}}{\int }{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{S}}_{21}\mathit{Cos}\mathit{\alpha }+\underset{\mathbf{B}}{\int }{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{S}}_{22}\mathit{Sin}\mathit{\alpha }+{\Omega }_1\hfill \\ \underset{\mathbf{B}}{\int }{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{S}}_{22}\mathit{Cos}\mathit{\alpha }-{\Omega }_1\underset{\mathbf{B}}{\int }{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{S}}_{21}\mathit{Sin}\mathit{\alpha }\}+{\mathit{k}}_2\mathit{Cm}\hfill \\ +{\mathit{k}}_2\mathit{Cm}\int [{\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{xS}}_{22}+{\mathit{C}}_{\mathit{p}}(\mathit{x}-{\mathit{x}}_{\mathit{ref}})\mathit{\delta }{\mathit{S}}_{22}-\hfill \\ {\mathit{C}}_{\mathit{p}}\mathit{\delta }{\mathit{yS}}_{2I}-{\mathit{C}}_{\mathit{p}}\mathit{y}\mathit{\delta }{\mathit{S}}_{2I}]{\Omega }_2\mathit{Cm}\underset{\mathbf{B}}{\int }[{\mathit{C}}_{\mathit{p}}\mathit{Sin}\mathit{\alpha }\mathit{\delta }{\mathit{S}}_{21}-\hfill \\ {\mathit{C}}_{\mathit{p}}\mathit{Cos}\mathit{\alpha }\mathit{\delta }{\mathit{S}}_{22}]+\underset{\mathbf{B}}{\int }{\mathit{C}}_{\mathit{p}}(\mathit{\delta }{\mathit{S}}_{21}{\mathit{\psi }}_2+\mathit{\delta }{\mathit{S}}_{22}{\mathit{\psi }}_3)\mathit{d}\mathbf{B}-\hfill \\ \underset{D}{\int }\frac{\mathit{\partial }{\mathit{\psi }}^{\mathit{T}}}{\mathit{\partial }{\mathit{\xi }}_{\mathit{k}}}\mathit{\delta }{\mathit{S}}_{\mathit{kl}}{\mathit{f}}_{\mathit{l}}\mathit{d}{D}_{\mathit{\xi }}\hfill \end{array}$

And the gradients can then be defined with respect to the design variables xi as:

(22) $\mathit{G}({\mathit{x}}_{\mathit{i}})=\frac{\mathit{\delta }\mathit{I}}{{\mathit{\delta }}_{{\mathit{X}}_{\mathit{i}}}}$

NUMERICAL METHOD FOR ADJOINT EQUATIONS

The semi discretized adjoint equations of Eq. 17 is obtained using the same approach used for the flow equations:

(23) $\begin{array}{c}{\forall }_{\mathit{ij}}\frac{\mathit{\partial }{\psi }_{\mathit{ij}}}{\mathit{\partial }\mathit{t}}-{\mathit{R}}_{\mathit{ij}}=0\hfill \\ {\mathit{R}}_{\mathit{ij}}=\frac{1}{2}\{{\mathit{C}}_{\mathit{l}\mathit{ij}}^{\mathit{T}}({\mathit{\psi }}_{\mathit{i}+1,\mathit{j}}-{\mathit{\psi }}_{\mathit{i}-1,\mathit{j}})+{\mathit{C}}_{2\mathit{ijk}}^{\mathit{T}}\hfill \\ ({\mathit{\psi }}_{\mathit{i},\mathit{j}+1}-{\mathit{\psi }}_{\mathit{i},\mathit{j}-1}\mathrm{)\}}-\{{\mathit{d}}_{\mathit{i}+\frac{1}{2},\mathit{j}}-{\mathit{d}}_{\mathit{i}-\frac{1}{2},\mathit{j}}+{\mathit{d}}_{\mathit{i},\mathit{j}+1/2}-{\mathit{d}}_{\mathit{i},\mathit{j}-1/2}\}\hfill \end{array}$

Introducing the artificial dissipation term to the Eq. 23 is crucially important to avoid discontinuities and to keep vector variable ψ differentiable. The five stage modified Runge-Kutta time stepping scheme is used to march the adjoint equations to the steady-state limit. Further details of the procedure of obtaining Eqs. 17 to 23 can be found in Jameson (1990, 1995).

OPTIMIZATION ALGORITHM

A simple optimization algorithm can be established by using a line search method. These methods require an algorithm to choose a direction p and search along this direction from the current design variables to obtain a new set of variables for the cost function value. Once the direction is chosen, then a step length λ is multiplied to the search direction to advance the optimization to the next iteration. A simple optimization algorithm can be defined by setting the search direction l=-G, to the negative of the gradient at every iteration:

(24) $\mathit{l}=-\mathit{G}$

With a line search method, the step size λ is chosen so that the maximum reduction of the objective function I(X) is attained. The search procedure used in this work is a descent method, in which small steps are taken in a direction defined by the smoothed gradient. X represents the design variable, and G=∇I the gradient. Instead of making the step:

(25) $\mathit{\delta }\mathit{X}=\mathit{\lambda }\mathit{l}=-\mathit{\lambda }G$

We replace the gradient ∇I by a smoothed value ∇I to guarantee the generation of a sequence of smooth shapes. When smoothing is performed in the x direction, the smoothed gradient ∇I may be calculated from a discrete approximation to:

(26) $\stackrel{─}{\nabla \mathit{I}}-\frac{\partial }{\partial \mathit{x}}\epsilon \frac{\partial }{\partial \mathit{x}}\stackrel{─}{\nabla \mathit{I}}=\nabla \mathit{I}$

where ε is the smoothing parameter.

Finding the new position of the boundary mesh points is also necessary, in order to estimate the variation of internal nodes accordingly. Jameson (1986, 1990) introduced a grid perturbation method that modifies the current location of the grid points based on perturbations at the geometry surface. This allows the variation of the grid point location in the equation for gradient evaluation, to be substituted with the variation of the surface points. The Optimization procedure can finally be summarized as follows:

1. Solve the flow equations for ρ, u ₁, u ₂, p and e.

2. Solve the adjoint equations for ψ subject to appropriate boundary conditions.

3. Evaluate the gradients $G=\frac{\mathit{\delta }\mathit{I}}{\mathit{\delta }{\mathit{X}}_{\mathit{S}}}$ and calculate the corresponding smoothed gradient G;

4. Update the shape based on the direction of steepest descent δX_s=- λ∇I;

5. Solve the adjoint equations again for determining the gradients $G=\frac{\mathit{\delta }\mathit{Cl}}{\mathit{\delta }{\mathit{X}}_{\mathit{S}}}$ and consequently for updating angle of attack by using the equation:

6. $\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathit{\alpha }}{|}_X\mathit{\delta }\mathit{\alpha }+\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathit{X}}{|}_{\mathit{a}}\mathit{\delta }{\mathit{X}}_{\mathit{S}}=0$

   or:

   (27) $\delta \mathit{\alpha }=-\frac{\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathit{X}}{|}_{\mathit{\alpha }}}{\frac{\mathit{\partial }\mathit{Cl}}{\mathit{\partial }\mathit{\alpha }}{|}_{\mathit{X}}}\mathit{\delta }{\mathit{X}}_{\mathit{S}}$

   where δX_s is the perturbed coordinates of nodes on the airfoil surface, obtained in step (d); and

7. Return to step (a) until convergence is reached.

RESULTS

In this section, we first validate our flow solver code by simulating several benchmark problems for the transonic inviscid flow over airfoils. For this purpose, transonic flow over NACA 0012 and RAE 2822 airfoils are considered. The validation is performed by comparing the obtained results to the previously published results by Pulliam et al. (1983) and Cook et al. (1979). In all cases, the far boundaries are located 5 chords away from the center of the airfoil. The NACA airfoil is set to three flow conditions: I) M_∞=0.63, α=2.0°; II) M_∞=0.8, α=0.0°; III) M_∞=0.8, α=1.25° and the free stream Mach number and the angle of attack for the RAE airfoil is IV)M_∞=0.67, α=1.06°.

The Courant number number was taken as 4.0 and five orders of magnitude reduction in the density residual is considered as the convergence criteria. In all cases, the convergence criterion is met in less than 1,500 time steps, as shown in Fig. 1.

Figure 1
Convergence History of solutions.

The computational grid consists of 128 × 33 nodes stretched near the leading and trailing edges. The grid is also clustered towards the airfoil surface. Computations were performed on a personal computer with a dual core CPU for which the computational time for the aforementioned cases is less than 2 minutes.

The pressure coefficients C_p on the upper and lower surface of the airfoils are shown in Figs. 2 to 5 and compared with the other computational results. In these Figures, X has been normalized and denotes the distance from the airfoil centre. The agreement between the Euler calculation and the other reference data is obvious. Figures 6 to 8 display the pressure contour plots for case III and IV. The smoothness of rapid expansion near the nose region on the upper surface is notable. The results prove that the region near stagnation point is calculated without any excessive numerical difficulty.

Figure 2
Comparison of pressure distribution for the case I.

Figure 3
Comparison of pressure distribution for the case II.

Figure 4
Comparison of pressure distribution for the case III.

Figure 5
Comparison of pressure distribution for the case IV.

Figure 6
Pressure contour plots for the case III.

Figure 7
Pressure contour plots for the case IV.

Figure 8
Inviscid drag minimization of NACA 0012, N =0.75, Cl =0.656, Initial Solution, Cd =0.0345, Cm =-0.0335, α =3.00º.

This problem is very important for the shape optimization procedure and can deteriorate the smoothness of surface. In all cases, the shock location is also computed accurately with a resolution of 2 to 3 grid points. Case III is a more difficult one, where a strong shock forms on the upper surface and a weak shock on the lower one. The grid is clustered near the shock position on the upper-lower surfaces of the airfoil for this case.

The performance of the optimization algorithm for the cost function is examined and the effects of weight coefficients k₁ and k₂ are explored. As our first test case, we set k₁=1.0 and k₂=0.0, therefore, reducing the problem to a pressure drag minimization one. NACA 0012 airfoil is considered to be the initial geometry and the flow variables are set to M_∞=0.75 and angle of attack α=3.0°.

The lift and drag coefficients for the initial geometry are calculated as Cl=0.656 and Cd=0.0345. For the initial geometry, a strong shock is developed on the upper surface of the airfoil. The vertical coordinates of the boundary grid points are considered as the design variables and the aim of optimization process is to reach a new airfoil geometry for which the drag coefficient is minimum while the lift coefficient is maintained at the same value of Cl=0.656. An implicit smoothing technique is also used to smooth out the gradients before modifying the location of each point on the airfoil surface.

The convergence is quickly attained by this optimizing technique and only after four design cycles the drag coefficient is reduced to one third of its initial value under the fixed lift coefficient constraint Cl = 0.656 . Comparisons of the pressure distribution for the new geometry in this figure and for the initial NACA 0012 indicate that the upper surface nodes are modified in a way that the local Mach number and the strength of the shock are attenuated on this surface. This, in turn, may lead to a lower wave drag but decreases the lift coefficient as well. As a result, the fall in the lift coefficient is compensated by modifying the shape of the trailing edge, as shown in Fig. 9 to recover the pressure difference near this region.

Figure 9
Inviscid drag minimization of NACA 0012, M=0.75, After 50 design iterations, Cl=0.656, Cd=0.0448, α =1.45º.

Figure 9 illustrates that the final design shaped is achieved after 50 designed cycles. The upper surface shock disappeared from the pressure distribution results and the drag coefficient is reduced to Cd=0.0448 , i.e., one eighth of its initial value. As can be seen in this Figure, the optimal shape is no longer symmetric, and a cusped shape is developed at the trailing edge.

The convergence history of the gradient norm and the drag coefficient are shown in Figs. 10 to 11. These Figures indicate that the convergence rate of the procedure is fast and that the optimization is accomplished mostly in the early cycles.

It is also required to show that the final result depends on the initial geometry. In order to do so, we repeat the computations for a five digit NACA 23012 airfoil as the initial geometry and change the angle of attack to reach the same lift coefficient as in the previous example.

To achieve this goal, the angle of attack must be changed to α=1.24°, as shown in Fig. 12. The initial drag coefficient is also calculated as Cd = 0.0433. Figure 13 illustrates that the optimization procedure converges to different solution for the final geometry and the pressure distributions and, therefore, the optimizing procedure is dependent to the initial conditions. We also notice that in this drag minimization problem, the pitching moment of the optimal shape is unfavorably increased from Cm=-0.0558 to Cm=0.154.

Figure 10
Convergence history of drag minimization of NACA 0012 airfoil.

Figure 11
Decrement of gradient norm of drag minimization of NACA 0012 airfoil.

Figure 12
Inviscid drag minimization of NACA 23012 M=0.75, Cl=0.660, Initial solution, Cd=0.0433, Cm=-0.0558, α =1.24º.

Figure 13
Inviscid drag minimization of NACA 23012, M=0.75, Final Solution, Cl=0.660 Cd=0.0436, α =-0.17º.

As our second example, the other extreme case of k₁=0.0 and k₂=1.0 is considered, i.e., the pitching moment coefficient is minimized under the constant lift constraint and without restricting the drag coefficient. The initial geometry and flow conditions are as in the previous example. As shown in Fig. 14, there is a 94% reduction in the pitching moment for the optimal shape, while the lift coefficient remains almost constant. The drag coefficient has also shown a 23% increase.

Figure 14
Pitching moment minimization of NACA 0012, M=0.75, After 36 design iterations, Cl=0.656, Cd=0.0424, Cm=-0.0002, α =3.71º.

To minimize the drag coefficient under constant lift and pitching moment coefficient, the weight constants are set to k₁=0.50 and k₂=0.50.

Figure 15 indicates the obtained results for the same initial airfoil geometry and the flow conditions. The obtained results revealed that the drag coefficient is reduced by 76% while the lift and pitching moment coefficients are remained constant. This is the optimal shape for an airfoil that can reduce the compressible pressure drag at a given required lift without compromising the pitching moment coefficient.

Figure 15
Optimization of NACA 0012, k₁ =0.50, k₂ =0.50,M=0.75, After 32 design iterations, Cl=0.656, Cd=0.0059, Cm=-0.0338, α =2.65º.

The test cases above show that the proposed cost function for the ASO is very efficient for designing airfoil with specific requirements in subsonic and transonic flows.

CONCLUSION

An efficient methodology for ASO of airfoils in transonic or high subsonic regime was described. In this approach, the Euler equations are solved for a symmetric and asymmetric airfoil with a multi-block method. The results obtained from this approach for the pressure distribution on the airfoil surface were validated and it was shown that the aerodynamic coefficients are estimated with a rather good accuracy. Then the prepared solver is employed to optimize the airfoils based on the adjoint method for minimum pressure drag, while maintaining the lift and pitching moment coefficients.

Then, a new cost function was introduced for the problem of reducing the pitching moment and drag coefficients under the constant lift constraint. The procedure was first tested for the extreme case of reducing pitching moment under the constant lift constraint without restricting the drag coefficient which was led to a final shape with higher wave drag coefficient. However, by proper selection of cost function weight coefficients, it is possible to reduce the wave drag under constant pitching moment and lift coefficients. In this case, the final solution is an airfoil with a sharp trailing edge and a narrow leading edge.

Drag minimization of airfoils, in the case of constant lift coefficient, has led to a cusped shape of trailing edge and relatively flat upper surface in both symmetric and cambered airfoils; but the final solution is different for two cases, especially for the shape of leading edge of the airfoils.

These results indicate that the optimization procedure based on the adjoint method depends on the initial shape of airfoil. In other words, this technique can converge to a local extremum and not a global one.

The main advantage of this method is that satisfactory results can be achieved by minimum computational efforts, particularly for the grid generation problem and the required central processing unit of computer (CPU) time. This feature becomes important where a low fidelity accurate case such as rotor blade flow solver is required for primarily design or optimization purposes. However, the proposed method ignores important phenomena such as the viscous flow or shock boundary layer interactions.

REFERENCES

• Anderson, W.K., Newman, J.C., Whitfield, D.L. and Nielsen, E.J., 2001, "Sensitivity analysis for Navier-Stokes equations on unstructured meshes using complex variables", AIAA Journal, Vol. 39, No. 1, pp. 56-63. doi:10.2514/2.1270.
  » https://doi.org/10.2514/2.1270
• Barger, R.L., 1975, "Procedures for the design of low-pitching-moment airfoils", NASA TN D-7982.
• Cook, P.H., Mc Donald, M.A. and Firmin, M.C.P., 1979, "Aerofoil RAE 2822 pressure distributions, and boundary layer and wake measurement", AGARD Advisory Report No. 138.
• Jameson, A., 1986, "Multigrid algorithms for compressible flow calculations", Proceedings of the 2nd European Conference on Multigrid Methods II, Lecture Notes in Mathematics, No.1228, pp. 166-201.
• Jameson, A., 1988, "Aerodynamic design via control theory", Journal of Scientific Computing, Vol. 3, No. 3, pp. 233-260.
• Jameson, A., 1990, "Automatic design of transonic airfoils to reduce the shock induced pressure drag", In Proceedings of the 31st Israel Annual Conference on Aviation and Aeronautics, pp. 5-17.
• Jameson, A., 1995, "Optimum aerodynamic design using CFD and control theory", AIAA Paper 95-1729-CP.
• Jameson, A. and Baker, T.J., 1983, "Solution of Euler equations for complex configurations", AIAA, pp. 83-1929.
• Jameson, A., Schmidt, W. and Turkel, E., 1981, "Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes", AIAA, pp. 81-1259.
• Hazra, S.B., Schulz, V., Brezillon, J. and , Gauger, N., 2005, "Aerodynamic shape optimization using simultaneous pseudo-timestepping", Journal of Computational Physics, Vol. 204, No. 1, pp. 46-64. doi:10.1016/j.jcp.2004.10.007.
  » https://doi.org/10.1016/j.jcp.2004.10.007
• Leifsson, L., Koziel, S. and Ogurtsov, S., 2011, "Inverse design of transonic airfoils using variable-resolution modeling and pressure distribution alignment", Procedia Computer Science, Vol. 4, pp. 1234-1243. doi:10.1016/j.procs.
  » https://doi.org/10.1016/j.procs
• 2011.04.133.Leifsson, L. and Koziel, S., 2010, "Multi-fidelity design optimization of transonic airfoils using shape-preserving response prediction", Procedia Computer Science, Vol.1, No. 1, pp. 1311-1320. doi:10.1016/j.procs.2010.04.146.
  » https://doi.org/10.1016/j.procs.2010.04.146
• Nadarajah, S. and Tatossian, C., 2008, "Multi-objective aerodynamic shape optimization for unsteady viscous flows", Journal of Optimization and Engineering, Vol. 11, No. 1, pp. 67-106. doi: 10.1007/s11081-008-9036-4.
  » https://doi.org/10.1007/s11081-008-9036-4
• Nadarajah, S. and Jameson, A., 2001, "Studies of the continuous and discrete adjoint approaches to viscous automatic aerodynamic shape optimization", AIAA 2001-2530, AIAA 15th. Computational Fluid Dynamic Conference, Anaheim, Ca, June 11-14.
• Pulliam, T.H., Jespersen, D.C. and Childs, R.E., 1983, "An enhanced version of an implicit code for the Euler equations", AIAA-83-0344. January 10-13, Reno, Nevada.
• Reuther, J.J., Jameson, A., Alonso, J.J., Rimlinger, M.J. and Saunders, D., 2001, "Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers", parts 1 and 2, Journal of Aircraft, Vol. 36, No. 1, pp. 51-74.
• Reuther, J. and Jameson, A., 1994, "Control based airfoil design using the Euler equations", AIAA paper 94-4272-CP.
• Xie, L., 2002, "Gradient-based optimum aerodynamic design using adjoint methods".
• Ying, G.Y., Feng, H. and Yu, S.M., 2011, "Aerodynamic airfoil design using the Euler equations based on the dynamic evolution method and the control theory", Science China Physics, Mechanics and Astronomy, Vol. 54, No. 4, pp. 697-702. doi:10.1007/s11433-011-4287-z.
  » https://doi.org/10.1007/s11433-011-4287-z
• Yu, L., Quan, C.X., Neng, L.Z. and Wu, X.J., 2010, "Multi-objective robust airfoil optimization based on non-uniform rational B-spline (NURBS) representation", Science China Physics, Technological Sciences, Vol. 53, No. 10, pp. 2708-2717. doi: 10.1007/s11431-010-4075-4.
  » https://doi.org/10.1007/s11431-010-4075-4

Publication Dates

• Publication in this collection
  Oct-Dec 2014

History

• Received
  04 Aug 2014
• Accepted
  14 Oct 2014