On an Exact Step Length in Gradient-Based Aerodynamic Shape Optimization

Abstract

This study proposeda novel exact expression for step length (size) in gradient-based aerodynamic shape optimization for an airfoil in steady inviscid transonic flows. The airfoil surfaces were parameterized using Bezier curves. The Bezier curve control points were considered as design variables and the finite-difference method was used to compute the gradient of the objective function (drag-to-lift ratio) with respect to the design variables. An exact explicit expression was derived for the step length in gradient-based shape optimization problems. It was shown that the derived step length was independent of the method used for calculating the gradient (adjoint method, finite-difference method, etc.). The obtained results reveal the accuracy of the derived step length.

1. Introduction

The combination of computational fluid dynamics with numerical optimization algorithms has become a numerical tool used to predict the optimal shape for aerodynamic configurations, such as airfoils. During the last few decades, aerodynamic shape optimization has been an important area of research and great advances have been made due to the rapid progress in the development of high-performance computers and computational algorithms. The traditional “trial and error” approach used to design aerodynamic configurations is based on wind-tunnel testing. Nowadays, the high cost associated with wind-tunnel testing of real aerodynamic models is significantly reduced using numerical shape optimization. Instead of costly experiments of a whole range of aerodynamic configurations, only the optimal shape predicted bythe computational analysis is tested in the wind tunnel. When validated this way, the optimal design may be used in the development of the final aerodynamic configuration. To make this approach viable, shape optimization algorithms for aerodynamic configurations must be as efficient as possible.

Among the factors involved in aerodynamic shape optimization are the type of flow regime (subsonic, transonic, supersonic, hypersonic), the governing equations (Laplace, Euler, Navier–Stokes), the grid used to mesh the flow domain (structured, unstructured, if structured: O-type, C-type, H-type, or hybrid), the design variables and the angle of attack, the method of discretization of the governing equations (finite element method (FEM), finite difference method (FDM), finite volume method (FVM)), the optimization method (Newton and quasi-Newton methods, conjugate gradient, genetic algorithm, etc.), the shape sensitivity analysis (finite-difference method, adjoint method, automatic differentiation, etc.), the aerodynamic body to be optimized(airfoil, wing, fuselage, blades, etc.), the objective function definition, the body parameterization method, and the imposed constraints (geometrical, etc.).

In Hicks and coworkers [1,2], computational fluid dynamics was combined with numerical optimization to predict theoptimal shape design for airfoils and wings. In Hicks et al. [1], the drag on nonlifting symmetric transonic airfoils in an inviscid flow is minimized under geometric constraints. The numerical optimization method is based on feasible directions and the inviscid flow is governed by the transonic, small-disturbance equation. In Hicks et al. [2], the optimal shape design of a wing is studied. The problems of drag minimization and lift-to-drag ratio maximization are addressed bysolving the three-dimensional full potential equation coupled with the conjugate gradient method. The first use of adjoint equations for optimal design problems in fluid dynamics is presented inPironneau [3] and for aerodynamic shape optimization problems governed by the full potential and the Euler equations arepresented in Jameson and coworkers [4,5,6,7]. The application of the adjoint method to aerodynamic shape optimization has been the subject of intensive research during the past few decades [8,9,10,11,12]. The optimal shape design in fluid mechanics and aerodynamics using different methods is thoroughly studied inMohammadi and Pironneau [13] and an excellent study on the coupling of computational fluid dynamics and optimization methods is given inThevenin and Janiga [14]. A thorough examination of various aerodynamic optimization methods is given inSkinner and Zare-Behtash [15]. An extensive study on the mathematical aspects of the discrete adjoint method is presented in Giles and coworkers [16,17]. InElliott and Peraire [18], the discrete adjoint method is used to calculate the gradient of the objective function with respect to design variables in airfoil and wing inverse design problems governed by the Euler equations using unstructured meshes. It is shown that the computation of the gradient of the objective function using the adjoint method is independent of thenumber of design variables. In Nemec and Zingg [19], a Newton–Krylov algorithm is applied to the discrete adjoint and the discrete flow-sensitivity approaches to calculate the gradient of the objective function in some two-dimensional aerodynamic shape optimization problems governed by the Navier–Stokes equations. The design problems considered include an inverse design, maximization of the lift-to-drag ratio, lift-constrained drag minimization, and lift enhancement. In Anderson and Venkatakrishnan [20], the continuous adjoint method is formulated and applied for aerodynamic shape optimization problems governed by the Euler and Navier–Stokes equations using unstructured grids for different objective functions, such as drag minimization, lift maximization, and an inverse design (based on a pressure distribution).

By increasing the number of design variables in large-scale aerodynamic shape optimization problems, the most efficient method forcalculating the involved gradients is to combine gradient-based optimization methods, such as steepest-descent, conjugate-gradient, or quasi-Newton method, with the adjoint method. In the gradient-based optimization methods, the iterative process used to update the solution is based on the calculation of two variables: (1) the direction of descent and (2) the step length (a positive scalar). The former, which may be obtained from the calculation of the gradient of the underlying objective function, specifies the direction along which the value of the objective function is reduced, while the latter determines the distance that the updated solution should move in the direction of descent to substantially minimize the objective function. The success of a gradient-based optimization method depends on effective choices of both the direction of descent and the step length [21]. In aerodynamic shape optimization problems, the gradient of the objective function may be determined usingthe finite-difference method or the adjoint method. The step length is calculated usingline search or trust-region methods. If the objective function to be minimized can be expressed as a polynomial function, one can find an exact step length; otherwise, an inexact step length may be obtained usingthe line search or trust-region methods, and it is chosen in such a way that some conditions, such as the Wolfe conditions, are satisfied [21]. To the best of our knowledge, in aerodynamic shape optimization problems, there exists no exact expression for the step length. In this study, an exact step length is developed for the first time for unconstrained minimization problems in aerodynamic shape optimization, such as drag and drag-to-lift ratio minimization problems.

2. Governing Equations

The steady inviscid transonic flows considered here are governed by the Euler equations. In this study, these equations were solved using the ANSYS Fluent v19.2 software (ANSYS, Inc., Canonsburg, PA, USA) to obtain the airfoil surface pressure distribution needed to calculate the sensitivity coefficients using the finitedifference method.

Airfoil Parameterization

In this study, Bezier curves (BezierAirfoil: A MATLAB code for generating an airfoil shape using Bezier curve-V2. DOI: 10.13140/RG.2.2.30002.56000) were used to parameterize the airfoil surfaces.

Bezier curve: The mathematical definition of a Bezier curve is as follows:

$$P(t)={\displaystyle \sum _{i=0}^n{B}_i{J}_{n,i}(t)}$$

(1)

where

$$J_{n,i}(t)=\frac{n!}{i!(n-i)!}{t}^i{(1-t)}^{n-i}$$

(2)

is the Bernstein basis polynomial of degree $n$, $t\in [0,1]$, and $P(t)=\left(x(t),y(t)\right)$. $x(t)$ and $y(t)$ are the $x$- and $y$- coordinates of the predetermined data points on the airfoil surface. By convention:$0^0\equiv 1$ and $0!\equiv 1$. The order of a Bezier curve is equal to $n+1$ (the number of the control points).

The parametric Bezier curve of interest here was of order 11 (because it can successfully generate a large number of airfoil shapes with a high degree of accuracy) and is given as follows:

$$\begin{array}{c}n=10\Rightarrow \mathrm{number} \mathrm{of} \mathrm{control} \mathrm{points}=11\\ P(t)={\displaystyle \sum _{i=0}^{10}{B}_i{J}_{10,i}(t)=}{B}_0{J}_{10,0}(t)+\dots +B_{10}{J}_{10,10}(t)\\ =B_0\frac{10!}{0!\left(10-0\right)!}{t}^0{(1-t)}^{10-0}+\dots +B_{10}\frac{10!}{10!\left(10-10\right)!}{t}^{10}{(1-t)}^{10-10}\end{array}$$

(3)

In this study, the $x$-coordinates of airfoil nodes were considered constant during the optimization, $P(t)=y(t)$. Therefore, for the predetermined $y$-coordinate of the airfoil nodes, we could write:

$$\begin{array}{c}{y}_{1,j}(t)=B_0{(1-t_{y_{1,j}})}^{10}+B_{1}10t_{y_{1,j}}{(1-t_{y_{1,j}})}^9+B_{2}45t_{y_{1,j}}^2{(1-t_{y_{1,j}})}^8+B_{3}120t_{y_{1,j}}^3{(1-t_{y_{1,j}})}^7\\ B_{4}210t_{y_{1,j}}^4{(1-t_{y_{1,j}})}^6+B_{5}252t_{y_{1,j}}^5{(1-t_{y_{1,j}})}^5+B_{6}210t_{y_{1,j}}^6{(1-t_{y_{1,j}})}^4+B_{7}120t_{y_{1,j}}^7{(1-t_{y_{1,j}})}^3\\ +B_{8}45t_{y_{1,j}}^8{(1-t_{y_{1,j}})}^2+B_{9}10t_{y_{1,j}}^9{(1-t_{y_{1,j}})}^1+B_{10}{t}_{y_{1,j}}^{10}\end{array}$$

(4)

Therefore, the Bezier control points $B_0,\dots ,B_{10}$ could be obtained at the first iteration for the predetermined $y$-coordinate of the initial airfoil. Two Bezier curves were needed to construct an airfoil surface: one for the upper surface (11 Bezier control points) and one for the lower surface (11 Bezier control points). Hence, we needed $2\times 11=22$ Bezier control points. $y_{1,j}$ are the $y$-coordinates of the airfoil nodes in which $j=1,\dots ,\frac{N+1}{2}$($N$ is the number of nodes on the airfoil surface) for the upper surface, $j=\frac{N+1}{2},\dots ,N$ for the lower surface (Figure 1), and $t_{y_{1,j}}$ is the value of $t\in [0,1]$ associated with the $y_{1,j}$ (Figure 2). The airfoil chord was divided into $\frac{N+1}{2}-1$ intervals of equal length $l$.

3. Aerodynamic Shape Optimization

3.1. Design Variables

Here, in the aerodynamic shape optimization, the design variables were the Bezier control points $B_{D{V}_i}$ ($i=1,\dots ,16$), as shown in Figure 3. The first and last two Bezier control points on the upper and the lower surfaces of the airfoil (hollow circles in Figure 3) were fixed and not considered to bedesign variables.For the ease of coding, the design variables were numbered consistent with the node numbering on the airfoil surface (see Appendix C). In other words, the design variables were numbered from right to left on the upper surface and from left to right on the lower surface (see Figure 3). The reason for fixing two control points adjacent to the airfoil trailing edge ($B_9^U$ and $B_9^L$) was that the values of these two control points were very small ($B_9^U=-1.35283\times {10}^{-5}$ and $B_9^L=1.35283\times {10}^{-5}$ for the NACA0012 airfoil, as shown in Figure 4), especially given that they should be multiplied by another small variable (finite-difference step-size: see Equation (17)), and even the double-precision calculations of the flow solver (here ANSYS Fluent) may not be able to detect the change of the pressure of the perturbed nodes. Thus, the computation of the sensitivity coefficients using the finitedifference method for these two control points would be problematic and would result in a zero sensitivity coefficient due to the subtraction of nearly equal terms. The complex-step derivative approximation is a remedy for this issue if one has access to the solver source code to substitute all real type variable declarations with complex declarations and define all functions and operators that are not defined for complex arguments [22]. However, in this study, the two control points mentioned above were fixed and are not considered as design variables.

3.2. Objective Function

The aerodynamic shape optimization problem of interest here consistedof the minimization of the drag ($D$) to lift ($L$) ratio, $\frac{D}{L}$, at a fixedangle of attack $\alpha$. This was an unconstrained optimization problem since no constraint was considered for it.

The objective function for this case is:

$$\mathcal{J}={(\frac{D}{L})}^2$$

(5)

where

$$D=A\cos \alpha +N\sin \alpha$$

(6)

$$L=-A\sin \alpha +N\cos \alpha$$

(7)

and $A$ and $N$ are axial and normal forces, respectively (Figure 5) [23]. We can write the drag and lift forces as follows:

$$D={\displaystyle \sum _{j=1}^{N-1}\left(p_{1,j}(y_{1,j}-y_{1,j+1})\cos \alpha +p_{1,j}(x_{1,j+1}-x_{1,j})\sin \alpha \right)}$$

(8)

$$L={\displaystyle \sum _{j=1}^{N-1}\left(-p_{1,j}(y_{1,j}-y_{1,j+1})\sin \alpha +p_{1,j}(x_{1,j+1}-x_{1,j})\cos \alpha \right)}$$

(9)

Therefore, the objective function becomes:

$$\min \mathcal{J}={(\frac{D}{L})}^2={(\frac{A\cos \alpha +N\sin \alpha }{-A\sin \alpha +N\cos \alpha })}^2={\left(\frac{{\displaystyle \sum _{j=1}^{N-1}\left(p_{1,j}(y_{1,j}-y_{1,j+1})\cos \alpha +p_{1,j}(x_{1,j+1}-x_{1,j})\sin \alpha \right)}}{{\displaystyle \sum _{j=1}^{N-1}\left(-p_{1,j}(y_{1,j}-y_{1,j+1})\sin \alpha +p_{1,j}(x_{1,j+1}-x_{1,j})\cos \alpha \right)}}\right)}^2$$

(10)

3.3. Sensitivity Analysis

The sensitivity analysis in gradient-based optimization is concerned with the computation of the derivative of the objective function with respect to the unknown variables. Suppose we wish to calculate the sensitivity of the objective function $\mathcal{J}$ defined by Equation (10) to the Bezier control points $B_{D{V}_i}$ ($i=1,\dots ,16$). This can be mathematically expressed as:

$$\frac{\partial \mathcal{J}}{\partial B_{D{V}_i}}=2(\frac{D}{L})\frac{\partial }{\partial B_{D{V}_i}}(\frac{D}{L})=2(\frac{D}{L})\frac{\frac{\partial D}{\partial B_{D{V}_i}}L-\frac{\partial L}{\partial B_{D{V}_i}}D}{L^2}=2(\frac{D}{L^3})(\frac{\partial D}{\partial B_{D{V}_i}}L-\frac{\partial L}{\partial B_{D{V}_i}}D)$$

(11)

The derivatives on the right-hand side of Equation (11) can be evaluated by taking the derivative of drag $D$, Equation (8), and lift $L$, Equation (9) with respect to the design variables $B_{D{V}_i}$, as follows:

$$\begin{array}{c}\frac{\partial D}{\partial B_{D{V}_i}}={\displaystyle \sum _{j=1}^{N-1}[\left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}(y_{1,j}-y_{1,j+1})+p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\right)\cos \alpha }+\\ \left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}(x_{1,j+1}-x_{1,j})+p_{1,j}\underset{=0}{\underbrace{\frac{\partial (x_{1,j+1}-x_{1,j})}{\partial B_{D{V}_i}}}}\right)\sin \alpha ]\end{array}$$

(12)

Factoring $\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}$, we obtain:

$$\frac{\partial D}{\partial B_{D{V}_i}}={\displaystyle \sum _{j=1}^{N-1}\left[\left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}\left((y_{1,j}-y_{1,j+1})\cos \alpha +(x_{1,j+1}-x_{1,j})\sin \alpha \right)\right)+p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\cos \alpha \right]}$$

(13)

In a similar fashion, we obtain:

$$\begin{array}{c}\frac{\partial L}{\partial B_{D{V}_i}}={\displaystyle \sum _{j=1}^{N-1}[\left(-\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}(y_{1,j}-y_{1,j+1})-p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\right)\sin \alpha }+\\ \left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}(x_{1,j+1}-x_{1,j})+p_{1,j}\underset{=0}{\underbrace{\frac{\partial (x_{1,j+1}-x_{1,j})}{\partial B_{D{V}_i}}}}\right)\cos \alpha ]\end{array}$$

(14)

Factoring $\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}$, we obtain:

$$\frac{\partial L}{\partial B_{D{V}_i}}={\displaystyle \sum _{j=1}^{N-1}\left[\left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}\left(-(y_{1,j}-y_{1,j+1})\sin \alpha +(x_{1,j+1}-x_{1,j})\cos \alpha \right)\right)-p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\sin \alpha \right]}$$

(15)

Thus, the derivative of the objective function, Equation (10) becomes:

$$\begin{array}{c}\frac{\partial \mathcal{J}}{\partial B_{D{V}_i}}=2(\frac{D}{L^3})(\frac{\partial D}{\partial B_{D{V}_i}}L-\frac{\partial L}{\partial B_{D{V}_i}}D)=\\ 2(\frac{D}{L^3})[{\displaystyle \sum _{j=1}^{N-1}\left[\left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}\left((y_{1,j}-y_{1,j+1})\cos \alpha +(x_{1,j+1}-x_{1,j})\sin \alpha \right)\right)+p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\cos \alpha \right]}L-\\ {\displaystyle \sum _{j=1}^{N-1}\left[\left(\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}\left(-(y_{1,j}-y_{1,j+1})\sin \alpha +(x_{1,j+1}-x_{1,j})\cos \alpha \right)\right)-p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}\sin \alpha \right]}D]=\\ 2(\frac{D}{L^3}){\displaystyle \sum _{j=1}^{N-1}[\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}}(y_{1,j}-y_{1,j+1})(L\cos \alpha +D\sin \alpha )+\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}(x_{1,j+1}-x_{1,j})(L\sin \alpha -D\cos \alpha )+\\ p_{1,j}\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}(L\cos \alpha +D\sin \alpha )]\end{array}$$

(16)

The sensitivity coefficients $\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}$ may be computed using the finitedifference method, which is chosen for the simplicity of implementation. Suppose we wish to calculate the sensitivity of the pressure of nodes on the airfoil surface, $p_{1,j}$$(j=1,\dots ,N-1)$, to the $y$-coordinate of the 16 Bezier control points ($B_{D{V}_i}$ ($i=1,\dots ,16$), 8 for each of the upper and lower surfaces)in this study. The sensitivity analysis can be performed by introducing a small perturbation to the $y$-coordinate of each of the Bezier control points. The grid generation (using perturbed values for the Bezier control points) and flow problem may be solved to obtain the new values for the pressure $p_{1,j}$. Using these values for the pressure, the dependency of the pressure $p_{1,j}$ on the perturbation of the $y$-position of Bezier control points can be evaluated. The central finitedifference method is used to compute these sensitivities as follows:

$$\frac{\partial p_{1,j}}{\partial B_{D{V}_i}}=\frac{p_{1,j}(B_{D{V}_i}+\epsilon B_{D{V}_i})-p_{1,j}(B_{D{V}_i}-\epsilon B_{D{V}_i})}{2\epsilon B_{D{V}_i}}$$

(17)

In this study, the finitedifference step-size was $\epsilon =0.005$. The term $\epsilon B_{D{V}_i}$ is the perturbation in the $y$-position of the Bezier control point $B_{D{V}_i}$. Since the sensitivity of each pressure $p_{1,j}$$(j=1,\dots ,N-1)$ to each $y$-position of the Bezier control point $B_{D{V}_i}$ ($i=1,\dots ,16$) is required, the computation of the sensitivity coefficients using this method requires 2 × (total number of control points − total number of fixed control points) = 2 × [2(11) − 2(3)] = 32 additional solutions of the flow problem at each iteration.

Moreover, the terms $\frac{\partial (y_{1,j}-y_{1,j+1})}{\partial B_{D{V}_i}}=\frac{\partial y_{1,j}}{\partial B_{D{V}_i}}-\frac{\partial y_{1,j+1}}{\partial B_{D{V}_i}}$ can be calculated using Equation (4). By having the value of the pressure at each node from the solution of flow equations (here, the Euler equations are used), the drag $D$ and the lift $L$ can be calculated using Equations (8) and (9), respectively. The angle of attack $\alpha$ is given a priori and is constant. The airfoil surface coordinates, $x_{1,j}$ and $y_{1,j}$, are also known from the Bezier curve equation (here, the $x_{1,j}$ are constant and do not change in the optimization process).

4. Optimization Method

In this study, the unconstrained nonlinear minimization problem was solved using Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (a quasi-Newton method). The quasi-Newton method is a gradient-based optimization method thatuses an approximation of the Hessian matrix or the inverse of a Hessian matrix [24]. In the BFGS method, an approximate inverse of the Hessian matrix, ${\mathcal{B}}^{\left(\mathrm{k}\right)}$, is updated iteratively until an optimal solution is found. The steps of the aerodynamic shape optimization using the BFGS method as implemented in this study are as follows.

(1)

Specify the physical domain; the boundary conditions; the problem conditions, such as the free stream Mach number; and the angle of attack.

(2)

Generate structured or unstructured grids over the domain. The airfoil surface is generated using aBezier curve of degree $n$.

(3)

Solve the flow problem regarding finding the pressure values at thegrid points of the physical domain and hence the airfoil surface.

(4)

Using Equations (8) and (9), compute the objective function ${\mathcal{J}}^{\left(\mathrm{k}\right)}$ Equation (10).

(5)

Compute the sensitivity coefficients using Equation (17).

(6)

Compute the gradient direction $\nabla {\mathcal{J}}^{\left(\mathrm{k}\right)}$ using Equation (16).

(7)

The initial approximate inverse of Hessian matrix, ${\mathcal{B}}^{(1)}$, is taken as the identity matrix, namely ${\mathcal{B}}^{(1)}=\mathbf{I}$.

(8)

Compute the direction of descent ${\mathbf{S}}^{\left(\mathrm{k}\right)}$ using the following equation:

$${\mathbf{S}}^{\left(\mathrm{k}\right)}=-{\mathcal{B}}^{\left(\mathrm{k}\right)}\nabla {\mathcal{J}}^{\left(\mathrm{k}\right)}$$

(18)

(9)

Compute the search step length ${\beta }^{\left(\mathrm{k}\right)}$. An explicit expression for computing the search step length will be presented later.

(10)

Evaluate the new values of the design variables (Bezier curve control points) using the following equation:

$${\mathbf{B}}_{DV}^{(\mathrm{k}+1)}={\mathbf{B}}_{DV}^{(\mathrm{k})}+{\beta }^{\left(\mathrm{k}\right)}{\mathbf{S}}^{\left(\mathrm{k}\right)}={\mathbf{B}}_{DV}^{(\mathrm{k})}-{\beta }^{\left(\mathrm{k}\right)}{\mathcal{B}}^{\left(\mathrm{k}\right)}\nabla {\mathcal{J}}^{\left(\mathrm{k}\right)}$$

(19)

(11)

Set the next iteration $(\mathrm{k}=\mathrm{k}+1)$ and return to step 2. Update the approximate inverse of the Hessian matrix as [24]:

$${\mathcal{B}}^{(\mathrm{k}+1)}={\mathcal{B}}^{\left(\mathrm{k}\right)}+(1+\frac{{({\mathbf{g}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathcal{B}}^{\left(\mathrm{k}\right)}{\mathbf{g}}^{\left(\mathrm{k}\right)}}{{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathbf{g}}^{\left(\mathrm{k}\right)}})\frac{{\mathbf{d}}^{\left(\mathrm{k}\right)}{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}}{{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathbf{g}}^{\left(\mathrm{k}\right)}}-\frac{{\mathbf{d}}^{\left(\mathrm{k}\right)}{({\mathbf{g}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathcal{B}}^{\left(\mathrm{k}\right)}}{{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathbf{g}}^{\left(\mathrm{k}\right)}}-\frac{{\mathcal{B}}^{\left(\mathrm{k}\right)}{\mathbf{g}}^{\left(\mathrm{k}\right)}{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}}{{({\mathbf{d}}^{\left(\mathrm{k}\right)})}^{\mathrm{T}}{\mathbf{g}}^{\left(\mathrm{k}\right)}}$$

(20)

where:

$${\mathbf{d}}^{\left(\mathrm{k}\right)}={\mathbf{B}}_{DV}^{(\mathrm{k}+1)}-{\mathbf{B}}_{DV}^{(\mathrm{k})}={\beta }^{\left(\mathrm{k}\right)}{\mathbf{S}}^{\left(\mathrm{k}\right)}$$

(21)

$${\mathbf{g}}^{\left(\mathrm{k}\right)}=\nabla {\mathcal{J}}^{(\mathrm{k}+1)}-\nabla {\mathcal{J}}^{\left(\mathrm{k}\right)}$$

(22)

There exists no exact search step length for highly nonlinear aerodynamic shape optimization problems. Therefore, an inexact step length using line search methods, such as backtracking, is determined and it is chosen in such a way that some conditions, such as Wolfe conditions, are satisfied [13,21,25]. Here, we present a novel method to obtain an exact step length in aerodynamic shape optimization problems involving the unconstrained minimization of common objective function definitions, such as drag and the drag-to lift-ratio. The formulation given here is concerned with the Bezier curve only but it can be extended to other airfoil parameterization methods.

Exact Step Length

The exact step length can be obtained as:

$${\beta }_{\mathrm{exact}}^{(\mathrm{k})}=\underset{\beta \ge 0}{\arg \min }\mathcal{J}({\mathbf{B}}_{DV}^{(\mathrm{k}+1)})=\underset{\beta \ge 0}{\arg \min }\mathcal{J}\left({\mathbf{B}}_{DV}^{\left(\mathrm{k}\right)}-{\beta }^{(\mathrm{k})}{\mathcal{B}}^{(\mathrm{k})}\nabla \mathcal{J}({\mathbf{B}}_{DV}^{\left(\mathrm{k}\right)})\right)$$

(23)

For simplicity, let $B_{{}_{D{V}_i}}^{(\mathrm{k})}=a_i^{(\mathrm{k})}$ and ${\mathcal{B}}^{(\mathrm{k})}\nabla {\mathcal{J}}^{(\mathrm{k})}=b_i^{(\mathrm{k})}$ ($i=1,\dots ,2n-2$), and suppose, for generality, that the Bezier control points associated with the trailing and leading edges are considered fixed (they are not considered as design variables). Thus, we can write:

$$B_{D{V}_i}^{(\mathrm{k}+1)}=a_i^{(\mathrm{k})}-{\beta }^{(\mathrm{k})}{b}_i^{(\mathrm{k})}$$

(24)

where $a_i^{(\mathrm{k})}$ is divided into $a_i^{U^{(\mathrm{k})}}$($i=1,\dots ,n-1$) for the upper surface and $a_i^{L^{(\mathrm{k})}}$($i=n,\dots ,2n-2$) for the lower surface. To obtain the exact step length ${\beta }_{\mathrm{exact}}^{{}^{(\mathrm{k})}}$, we minimize $\mathcal{J}(B_{D{V}_i}^{(\mathrm{k}+1)})={\left(\frac{D(B_{D{V}_i}^{(\mathrm{k}+1)})}{L(B_{D{V}_i}^{(\mathrm{k}+1)})}\right)}^2$ with respect to ${\beta }^{(\mathrm{k})}$ as follows:

$$\frac{\partial \mathcal{J}}{\partial \beta }=2\frac{D}{L}\frac{\partial }{\partial \beta }(\frac{D}{L})=0\Rightarrow D=0$$

(25)

$$\begin{array}{c}D={\displaystyle \sum _{j=1}^{N-1}\left(y_{1,j}(p_{1,j}-p_{1,j-1})\cos \alpha +N\sin \alpha \right)}=\\ \left({\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}{y}_{1,j}(p_{1,j}-p_{1,j-1})+{\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}{y}_{1,j}(p_{1,j}-p_{1,j-1})}}\right)\cos \alpha +N\sin \alpha =0\end{array}$$

(26)

Since $y_{1,1}=y_{1,N}=y_{1,\frac{N+1}{2}}=0$, the indices in Equation (26) do not include $j=1,\frac{N+1}{2},N$. Equation (26) is the same as Equation (8) with the rearrangement of terms being with respect to $y_{1,j}$ rather than $p_{1,j}$. From the definition of the Bezier curve, Equations (1) and (2)), we can write the drag, Equation (26) as:

$$\begin{array}{c}D=\{{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}\left({\displaystyle \sum _{i=1}^{n-1}{B}_{{}_{D{V}_i}}^{U^{(\mathrm{k}+1)}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})+\\ {\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}\left({\displaystyle \sum _{i=1}^{n-1}{B}_{{}_{D{V}_i}}^{L^{(\mathrm{k}+1)}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})\}\cos \alpha +N\sin \alpha =0\end{array}$$

(27)

Using $B_{D{V}_i}^{(\mathrm{k}+1)}=a_i^{(\mathrm{k})}-{\beta }^{(\mathrm{k})}{b}_i^{(\mathrm{k})}$ (Equation (24)), we can write:

$$\begin{array}{c}D=\{{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}\left({\displaystyle \sum _{i=1}^{n-1}(a_i^{U^{(\mathrm{k})}}-{\beta }^{(\mathrm{k})}{b}_i^{U^{(\mathrm{k})}})\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})+\\ {\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}\left({\displaystyle \sum _{i=1}^{n-1}(a_i^{L^{(\mathrm{k})}}-{\beta }^{(\mathrm{k})}{b}_i^{L^{(\mathrm{k})}})\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})\}\cos \alpha +N\sin \alpha =0\end{array}$$

(28)

Afterexpanding, we obtain:

$$\begin{array}{c}D=\{{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}\left({\displaystyle \sum _{i=1}^{n-1}{a}_i^{U^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}-{\beta }^{(\mathrm{k})}{\displaystyle \sum _{i=1}^{n-1}{b}_i^{U^{(\mathrm{k})}}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})+\\ {\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}\left({\displaystyle \sum _{i=1}^{n-1}{a}_i^{L^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}-{\beta }^{(\mathrm{k})}{\displaystyle \sum _{i=1}^{n-1}{b}_i^{L^{(\mathrm{k})}}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})\}\cos \alpha +N\sin \alpha =0\end{array}$$

(29)

Using the identity ${\displaystyle \sum _{i=k_0}^{k_1}{\displaystyle \sum _{j=l_0}^{l_1}{a}_{ij}=}}{\displaystyle \sum _{j=l_0}^{l_1}{\displaystyle \sum _{i=k_0}^{k_1}{a}_{ij}}}$, we can write the above expression as:

$$\begin{array}{c}D=\{{\displaystyle \sum _{i=1}^{n-1}\left({\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}{a}_i^{U^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}-{\beta }^{(\mathrm{k})}{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}{b}_i^{U^{(\mathrm{k})}}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})+\\ {\displaystyle \sum _{i=1}^{n-1}\left({\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}{a}_i^{L^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}-{\beta }^{(\mathrm{k})}{\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}{b}_i^{L^{(\mathrm{k})}}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}}\right)}(p_{1,j}-p_{1,j-1})\}\cos \alpha +N\sin \alpha =0\end{array}$$

(30)

Therefore, we obtain

$${\beta }_{\mathrm{exact}}^{(\mathrm{k})}=\frac{\left[{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}{a}_i^{U^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}(p_{1,j}-p_{1,j-1})+{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}{a}_i^{L^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}(p_{1,j}-p_{1,j-1})}}}}\right]\cos \alpha +N\sin \alpha }{\left[{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=2}^{\frac{N+1}{2}-1}{b}_i^{U^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}(p_{1,j}-p_{1,j-1})+{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=\frac{N+1}{2}+1}^{N-1}{b}_i^{L^{(\mathrm{k})}}\frac{n!}{i!(n-i)!}{t}_{y_{1,j}}^i{(1-t_{y_{1,j}})}^{n-i}(p_{1,j}-p_{1,j-1})}}}}\right]\cos \alpha }$$

(31)

which is the exact step length. As can be seen from Equation (31), the value of the exact step length depends on the gradient of the objective function and it is independent of the method used for calculating the gradient. In other words, one can calculate the gradient with any available method, such as the adjoint method, the finitedifference method, etc., and then use Equation (31) to calculate the exact step length. It can be concluded that an inaccuracy in computing the gradient of the objective function or the approximate inverse of Hessian matrix ${\mathcal{B}}^{\left(\mathrm{k}\right)}$ can give rise to an inaccuracy in calculating the exact step length.

For the sake of simplicity, the derivation of the above general step length expression is presented for the Bezier curve of degree 10 and is given in Appendix A and a Fortran code to calculate it is given in Appendix D. The derivation of the exact step length for other forms of the objective function, such as drag minimization and lift-constrained drag minimization, is given in Appendix B.

5. Results

As mentioned previously, the objective of the aerodynamic shape optimization problem considered in this study was to minimize the drag-to-lift ratio (which is equivalent to the maximization of the lift-to-drag ratio) for an airfoil moving in a steady inviscid transonic flow through the implementation of the exact step length in a gradient-based optimization algorithm. The following test case is presented to investigate the accuracy and performance of the derived step length.

Test Case: A NACA 0012 airfoil operating at a fixed angle of attack of $2^{\circ }$ and a free-stream Mach number of 0.75 was considered for this test case as theinitial shape and flow conditions. The free stream pressure was $p_{\infty }=\mathrm{101,325.0} \mathrm{Pa}$. There was a strong shock at half of the airfoil chord length on the upper surface. Since there was shock-induced pressure drag on the airfoil inthese operating conditions, we aimedto optimize the airfoil to eliminate or weaken the shock as much as possible. The initial drag and lift forces were $D=515.760 \mathrm{N}$ and $L=\mathrm{17,090.0} \mathrm{N}$, respectively (this means that for the initial airfoil, NACA 0012, inthese flow conditions,$\frac{L}{D}=\frac{\mathrm{17,090.0} \mathrm{N}}{515.760 \mathrm{N}}=33.1$). The initial drag and lift coefficients were $c_d=1.294\times {10}^{-2}$ and $c_l=0.42867$, respectively, which are in excellent agreement with $c_d=1.25\times {10}^{-2}$ and $c_l=0.4225$ in Carpentieri et al. [26]. As shown in

Table 1. Results for the test case.

$\mathbf{Iteration} \mathbf{Number} (\mathbf{k})$	$\mathit{J}={(\frac{\mathit{D}}{\mathit{L}})}^2$	$\mathbf{Step} \mathbf{Length} \mathit{\beta }({\mathbf{m}}^2)$	$\mathit{D}(\mathbf{N})$	$\mathit{L}(\mathbf{N})$	${\mathit{c}}_{\mathit{d}}$	${\mathit{c}}_{\mathit{l}}$
0 (initial shape)	9.10775 × 10−4	-	515.760	17,090.000	1.294 × 10−2	0.42867
1	1.97841 × 10−5	0.5534	34.032	7651.200	8.536 × 10−4	0.19191
2	1.34569 × 10−5	1.6872	35.644	9716.600	8.940 × 10−4	0.24372
3	1.06817 × 10−5	1.1510	38.239	11,700.000	9.591 × 10−4	0.29348
4	9.92128 × 10−6	2.4779	40.664	12,910.000	1.020 × 10−3	0.32381
5	8.40409 × 10−6	0.0188	40.980	14,136.000	1.028 × 10−3	0.35456
6	4.89131 × 10−6	−1.2468	35.521	16,061.000	8.910 × 10−4	0.40284
7	1.50193 × 10−5	47.6144	72.886	18,807.000	1.828 × 10−3	0.47174
8	2.14748 × 10−5	2.8478	49.214	10,620.000	1.234 × 10−3	0.26638
9	1.66016 × 10−5	0.3450	48.401	11,879.000	1.214 × 10−3	0.29795
10	1.24568 × 10−5	0.2160	45.653	12,935.000	1.145 × 10−3	0.32445
11	7.33858 × 10−6	0.2599	38.305	14,140.000	9.608 × 10−4	0.35466
12	6.68139 × 10−6	0.0009	40.202	15,553.000	1.008 × 10−3	0.39011
13	8.65189 × 10−6	11.8508	40.162	13,654.000	1.007 × 10−3	0.34247
14	1.26244 × 10−5	1.8642	59.411	16,721.000	1.490 × 10−3	0.41941
15	5.73953 × 10−6	0.9977	36.142	15,086.000	9.065 × 10−4	0.37839
-	99.4% reduction	-	93.0% reduction	11.7% reduction	-	-

, the drag and lift forces of the optimal shape were $D=36.142 \mathrm{N}$ and $L=\mathrm{15,086.0} \mathrm{N}$, respectively (this means that for the optimal airfoil,$\frac{L}{D}=\frac{\mathrm{15,086.0} \mathrm{N}}{36.142 \mathrm{N}}=417.4$, which was an approximately 11.6 times increase in the value of ($\frac{L}{D}$). A substantial decrease in the drag (from $515.760 \mathrm{N}$ to $36.142 \mathrm{N}$) using only 16 design variables was achieved. Moreover, the drag and lift coefficients of the optimal shape were $c_d=9.065\times {10}^{-4}$ and $c_l=0.37839$, respectively. A structured O-grid (see Figure 6) with $N=601$ nodes on the airfoil surface and 109 nodes in the normal direction to the airfoil surface was generated. The inviscid transonic flow equations (Euler equations) were solved using the ANSYS Fluent software to compute the pressure at each node on the airfoil surface (for unperturbed and $2\times 16=32$ perturbed shapes foreach iteration). The computed pressures were then substituted into Equation (17) to calculate the sensitivity coefficients needed to compute the gradient of the objective function. A comparison of the initial shape (NACA 0012 airfoil) and the optimal one, as well as the corresponding surface pressure coefficients, is shown in Figure 7a,b, respectively. As can be seen, the shock on the initial shape was eliminated in the optimal shape. The static pressure distributions of the initial and optimal shapes are depicted in Figure 8a,b, respectively. A comparison of the initial and updated shape at each iteration is shown in Figure 9. As can be seen in Figure 9, the shock waseliminated or significantly weakened at each iteration (it is apparent from the value of the drag force given in Table 1, as well as from the thinner leading edge in the updated airfoils). The optimal shape is also in good agreement with the optimal shape in References [26,27]. A comparison of the initial shape (NACA 0012) and the optimal shapes obtained in Chen et al. [27] and from our study is shown in Figure 10. As can be seen, the upper surface shapes in both optimal shapes wereshock-free, which resulted in the minimum drag in both shapes (as can be seen in Table 1, the majority of the reduction in the objective function ${(\frac{D}{L})}^2$ in the test case was due to the decrease in the drag value). The history of the objective function is shown in Figure 11. The value of the objective function was substantially reduced by 99.4% and most of the reduction was obtained in the first iteration. The optimization process was continued for 15 iterations. The value of the exact step length used in the BFGS optimization method (using Equation (31) or Equation (A5)) is plotted in Figure 12. At iteration 6, the negative value for the step length was obtained (${\beta }^{(6)}=-1.2468$). The minimum value for the objective function was obtained at iteration 6; however, as this value was obtained by a negative step length, the optimization process was continued further. The sensitivity coefficients at some arbitrary nodal points (node numbers 20, 170, and 270 on the upper surface and node numbers 350, 500, and 580 on the lower surface; there are $\frac{N+1}{2}=\frac{601+1}{2}=301$ nodes on each surface of the airfoil) on the airfoil surface using three different finite-difference step-sizes $\epsilon ={10}^{-2}$, $\epsilon ={10}^{-3}$, and $\epsilon =5\times {10}^{-3}$ are depicted in Figure 13. As shown, the sensitivity coefficients at each node using the three step-sizes were in good agreement. However, in this study, the finite-difference step-size $\epsilon =5\times {10}^{-3}$ was used to calculate the sensitivity coefficients using Equation (17). The gradients of the objective function $\mathcal{J}={(\frac{D}{L})}^2$ with respect to the design variables $B_{D{V}_i}$ ($i=1,\dots ,16$) for the first and last iterations(iterations 1 and 15) are shown in Figure 14a,b, respectively, representing a significant decrease in the norm of the gradient (from $0.0697$ to $4.151\times {10}^{-4}$; a $99.4\%$ reduction) to infer the optimal shape.

It should be emphasized that the gradient of the objective function with respect to the design variables was computed using the finitedifference method. More accurate results are expected if the gradient computation is accomplished with more accurate methods, such as the adjoint method or the complex-step method.

6. Conclusions

This study aimed to develop a novel and exact step length for gradient-based aerodynamic shape optimization problems. The steady inviscid transonic flow governed by the Euler equations was considered and the Euler equations were solved usingthe finite element solver ANSYS Fluent. Minimization of the drag-to-lift ratio (which is equivalent to the maximization of the lift-to-drag ratio) was performed. The quasi-Newton (BFGS) method was employed to minimize the defined objective function. The airfoil (NACA 0012) surface was parameterized usinga Bezier curve and an explicit exact step length was derived for the unconstrained minimization problem. We showed that the value of the exact step length dependedon both the values of the gradient of the objective function and the Bezier control points. A test case was presented to investigate the effect of the exact step length on the shape optimization process. Using the exact step length, the optimal shape was obtained within a few iterations with a substantial reduction in the objective function value.