Abstract
This paper is concerned with an optimal shape design problem in aerodynamics. The inverse problem in question consists in finding the optimal shape an airfoil placed in a potential flow at a given angle of attack should have such that the pressure distribution on its surface matches a desired one. The numerical method to achieve this aim is based on a body-fitted grid generation technique (elliptic, O-type) to generate a mesh over the airfoil surface and solve for the flow equation. The O-type scheme is used due to its ability to generate a high quality (fine and orthogonal) grid around the airfoil surface. This paper describes a novel and very efficient sensitivity analysis scheme to compute the sensitivity of the pressure distribution to variation of grid node positions and both the conjugate gradient method (CGM) and a version of the quasi-Newton method (i.e., BFGS) are used as optimization algorithms to minimize the difference between the computed pressure distribution on the airfoil surface and desired one. The elliptic grid generation technique allows us to map the physical domain (body) onto a fixed computational domain and to discretize the flow equation using the finite difference method (FDM).
1. Introduction
Thanks to the advent of modern high speed computers over the last few decades, computational fluid dynamics (CFD) has been extensively employed as an analysis and as a design optimization tool. Among the methodologies often employed in shape optimization are gradient-based techniques. These techniques may be applied to minimize a specified objective function. In airfoil shape optimization, the objective function can be, for example, a measure of difference between the pressure distribution on the airfoil surface and a desired one, and it would be desirable to minimize this objective function. In this paper, we consider the 2D shape optimization of an airfoil in an irrotational and incompressible flow governed by the Laplace equation. The procedure employed is based on the elliptic grid generation, a novel sensitivity analysis (based on finite difference method), and an optimization method. The conjugate gradient method and an efficient version of quasi-Newton method (BFGS) will be used as the optimization algorithms. The airfoil surface is parameterized using the grid points and the Bezier curve. Three different types of design variables were considered: the grid points, the Bezier curve control points, and the maximum thickness of NACA00xx airfoils. It will be represented that the use of the Bezier curve significantly improves the optimization performance to reach the optimal shape. Furthermore, it will be shown that the proposed sensitivity analysis method reduces the computation cost significantly even for large number of the design variables.
Some of the earliest studies using a combination of CFD with numerical optimization in aerodynamic were made by Hicks et al. [1] and Hicks and Henne [2]. In [1], a procedure for optimal design of symmetric low-drag, nonlifting transonic airfoils in inviscid flow is proposed. The proposed procedure uses an optimization program based on the method of the feasible directions coupled with an analysis program that utilizes a relaxation method to solve the partial differential equation that governs the inviscid, transonic, and small disturbance fluid flow. The drag minimization with geometric constraints is considered in this reference. In fluid dynamics, Pironneau was the first one to use the adjoint equations for design [3]. This is the first application of control theory to design optimization. However, within the field of aeronautical computational fluid dynamics, Jameson was the first researcher who used the continuous adjoint formulation for aerodynamic shape optimization in transonic potential flows and flows governed by Euler equations [4–7]. Giles et al. made considerable contributions to the development of the discrete adjoint approach [8–11]. In [10], the adjoint equations are formulated for the transonic design applications for which there are shocks. The adjoint equations were already formulated for the incompressible or subsonic flows in which the assumption that the original nonlinear flow solution is smooth is valid. In [11], a number of algorithm developments are presented for adjoint methods using the “discrete” approach. In continuous adjoint method, the original partial differential equations are linearized, the adjoint partial differential equation and appropriate boundary conditions are formulated, and finally the equations are discretized. Unlike the continuous adjoint approach, in discrete adjoint approach the partial differential equations are discretized, the discrete equations are linearized, and then the transpose of the linear operator is used to form the adjoint problem. The adjoint equations have also been used by Baysal and Eleshaky to infer the optimal design for a scramjet-afterbody configuration which yields the maximum axial thrust [12] and by Ta’asan et al. to obtain an optimal airfoil shape [13]. Baysal and Eleshaky’s work was based on a computational fluid dynamics-sensitivity analysis algorithm (two different quasi-analytical approaches: the direct method and the adjoint variable method) to solve Euler equations for the inviscid analysis of the flow. Adjoint methods have been applied to incompressible viscous flow problems by Cabuk et al. [14] and Desai and Ito [15]. Cabuk et al. worked on the problem of determining the profile of a channel or duct that provides the maximum static pressure rise by solving the incompressible, laminar flow governed by the steady state Navier-Stokes equations. Early applications of discrete adjoint methods on unstructured meshes can be found in works by Elliott and Peraire in inviscid [16] and viscous flows [17] for 2D and 3D [18] configurations. In [16], an inverse design procedure for single- and multielement airfoils using unstructured grids and based on the Euler equations is presented. The discrete adjoint method is used to compute the sensitivities and the results are compared with corresponding finite difference values. It is shown that the use of the adjoint method practically eliminates the dependence of the objective function gradient computation on the number of design variables. The continuous adjoint approach for unstructured grids has been developed by Anderson and Venkatakrishnan [19]. In [19], aerodynamic shape optimization on unstructured grids using a continuous adjoint approach is developed and analyzed for inviscid and viscous flows. B spline and Bezier curves are employed to parameterize the airfoil surface. The objective functions considered include drag minimization, lift maximization, and matching a specified pressure distribution. The quasi-Newton optimization method is used to obtain the optimal design. Evolutionary algorithms, as methods that do not need the computation of the gradient, have recently gained much attention in the context of aerodynamic shape optimization [20–24]. Although they are of extremely high computational cost, they have the advantage that they can escape from a “local minimum” (a major issue in using gradient based methods) and have the ability of finding globally optimum solutions amongst many local optima [25, 26]. A detailed study of many methods in shape optimization in fluid mechanics is given by Mohammadi and Pironneau [27].
The adjoint approach, as an alternative to the finite difference method to compute the gradient of functional with respect to the design variables, is computationally very efficient. Therefore, as far as the computational cost is concerned, it is the appropriate choice. This is the case when there are a large number of design variables which makes use of the finite difference method impractical. The differences between the adjoint method and the finite difference one (to compute the gradient of functional with respect to the design variables) can be summarized as follows.
Adjoint Method. design variables, flow solution, and adjoint calculation.
Finite Difference Method. design variables, flow solutions. Because where is the objective function and are the design variables [28]. As can be seen, aerodynamic shape optimization with large number of design variables is computationally practical only when the adjoint method is used. However, as will be shown in this paper, a novel sensitivity analysis will be presented which makes use of the finite difference method comparable (from computation cost viewpoint) to the use of adjoint method. The numerical algorithm used in this paper is already employed in shape optimization problems in heat conduction [29, 30]. The numerical algorithm consists of three steps, namely, grid generation and flow equation solver to find the pressure on the airfoil surface, sensitivity analysis to compute the gradient of the objective function with respect to the design variables, and an optimization method to minimize the functional and reach optimum solution.
2. Governing Equation
For a two-dimensional incompressible flow, a stream function can be defined such that where and are the components of the velocity vector ; that is, ( and are the unit vectors in and directions, resp.). Combining the above definitions with the irrotationality condition leads to the following Laplace equation for the stream function
Consider an irrotational incompressible flow over an airfoil (Figure 1). The boundary conditions are as shown in Figure 1.
Conditions at Infinity. Far away from the airfoil surface (toward infinity), in all directions, the flow approaches the uniform freestream conditions. If the angle of attack (AOA) is , the free stream velocity , the components of the flow velocity can be written as
Condition on the Airfoil Surface. The relevant boundary condition at the airfoil surface for this inviscid flow is the no-penetration boundary condition. Thus the velocity vector must be tangent to the surface. This wall boundary condition can be expressed by where is tangent to the surface. In the problem of the flow over an airfoil, if the free stream velocity and the angle of attack are known, from the boundary conditions at infinity (see (4)) and the wall boundary condition (see (5)) one can compute the stream function at any point of the physical domain (flow region). Then, by knowing , one can compute the velocity of all points. Since for an incompressible flow, the pressure coefficient is a function of the velocity only, one can obtain the pressure of any point in the flow region, as will be shown.
Pressure Coefficient. The pressure coefficient is defined as where is the velocity of fluid at the point at which the pressure coefficient is being evaluated. At standard sea level conditions, where and are the freestream density and pressure, respectively. From (6),
indicates that the point at which the pressure coefficient is being evaluated is located at infinity.
indicates that the point at which the pressure coefficient is being evaluated is a stagnation point (where ). For an incompressible flow, this is the maximum allowable value of anywhere in the flow field.
And in regions of the flow where , value will be negative.
3. Grid Generation and Flow Solver
To calculate the pressure at any point in the flow region, a grid should be generated over the region. The grid generation method considered in this study is the elliptic grid generation, which was proposed by Thompson et al. [31] and is based on solving a system of elliptic partial differential equations to distribute nodes in the interior of the physical domain by mapping the irregular physical domain from the and physical plane (Figure 2) onto the and computational plane (Figure 3), which is a regular region.
The O-type elliptic grid generation technique is employed here which results in a smooth and orthogonal grid over the airfoil surface. The discretization of the physical domain (flow region) and the corresponding computational domain are shown in Figures 2 and 3, respectively. In the computational domain, and are the number of nodes in the and directions, respectively. The resulting O-type gird scheme over an airfoil for the case and or is shown in Figure 4.
The initial guess for the elliptic grid generation is performed using the transfinite interpolation (TFI) method. Since TFI method is an algebraic technique and does not require much computational time, it will be an appropriate initial guess for the elliptic grid generation method and accelerates convergence time for the elliptic method. Another advantage of using the TFI method as an initial guess is that it prevents the grids generated by the elliptic (O-type) method from folding.
If and are known, then from (4) one can obtain the stream function at any point on the boundaries of the physical domain as follows: where subscripts and refer to any two arbitrary grid points on the boundaries of the physical domain. Equations (8) and (9) are applied to vertical and horizontal boundaries of the physical domain, respectively. By knowing the values of the stream function on the boundaries of the physical domain as well as on the airfoil surface, we can obtain the values of over the physical domain by applying the Kutta condition [32, 33] and using the following formula (by mapping the physical domain onto the computational domain [29]): where and are grid control functions which control the density of grids towards a specified coordinate line or about a specific grid point. Equations (10) and (11) are discretized using the finite difference method. For more details, please refer to [29].
Velocity Calculation. There are three sections where the velocity must be known:(1)the outer boundaries (four sides , , , and of the rectangle shown in Figure 2);(2)the airfoil surface ( in Figure 2);(3)the inside of the physical domain.
The velocity values on the outer boundaries are known from the conditions at infinity (using (4)). In other words, -component of the velocity vector () on all the outer boundaries is equal to and -component of the velocity vector () on all the outer boundaries is equal to . For the inside of the physical domain and the airfoil surface, we can use the flowing relationships to evaluate the velocity. These relationships are obtained by using the transformation relationships and chain rule in mapping the physical domain onto the computational one. Consider The central and forward difference schemes are used for the inside of the physical domain and the airfoil surface, respectively. After obtaining the components of the velocity vector, the total velocity (velocity distribution) can be computed by As stated before, for an incompressible flow, the pressure coefficient can be expressed in terms of the velocity only. Thus (6) can be used to determine the pressure at any grid point in the domain. Therefore,
Validation of the Results for the Pressure Distribution. The results obtained here are compared with the results given in [34] which are obtained both analytically and by using the panel method (see Figure 7).
Validation Case. The pressure coefficient distribution () over the NACA 0012 airfoil at an angle of attack is plotted. The results are compared with the results from [34]. The O-type grid size used in the computation is 155 × 155. The computation time is 53 seconds.
4. Airfoil Parameterization
So far, the airfoil surface is parameterized by grid points which result in accurate pressure distribution on the airfoil surface (see Figures 14, 19, and 25). However, a large number of grid points are needed to obtain such accurate results which in turn lead to high (see Figures 5 and 6) computation cost. The design variables are the coordinate (usually -coordinate) of grid points. Therefore, the optimization process may be inappropriate if there are a large number of design variables since it is difficult to maintain a smooth geometry, the optimization problem will be difficult to solve, and the optimization strategy is likely to fail or be impractical [35]. Thus alternative methods of airfoil surface parameterization are needed. These methods should represent great flexibility in defining the airfoil surface with minimum design variables. In this paper, in addition to the grid points to represent the airfoil surface, Bezier curves (a special subset of B-spline) are employed due to their ability to produce airfoil surfaces easily and precisely with only a few control points.
Bezier Curve. A Bezier curve is a special case of a B-spline curve and is mathematically defined by where is Bernstein basis polynomial of degree . By convention and . Here, , the degree of the Bernstein basis polynomial is one less than the number of points in the Bezier polygon. In other words, the number of control points is . The points are the vertices of a Bezier polygon or the control points of a Bezier curve. The curve begins at and ends at . The order of a Bezier curve is equal to . In other words, the order of a Bezier curve is equal to the number of the control points [36].
In this paper, two different Bezier curves of order 7 (degree = 6) and of order 11 (degree = 10) will be considered. As it will be shown, the Bezier curve of order 7 represents the better optimization performance due to its less design variables. However, this kind of Bezier curve is not able to produce very accurate airfoil shapes. Indeed, it is appropriate to NACA 00xx airfoils only. On the other hand, the Bezier curve of order 11 can successfully generate any airfoil shape with a high degree of accuracy. Therefore, the formulation for the Bezier curve of order 11 only will be given here. The formulation for the Bezier curve of order 7 can be written in a similar fashion.
The parametric Bezier curve of order 11 is as follows: Therefore, In order to construct the airfoil surface, two Bezier curve will be considered corresponding to the upper and lower surfaces, respectively. Here there are 11 control points (vertices) for each surface. Since the coordinates of the airfoil surface are known, the problem is to determine values for the control points . In other words, our problem is to specify the coordinates of the control points so that the curve passes through the predetermined data points on the airfoil surface. Equation (16) can be written in matrix form as follows: If the number of the chosen data points on the airfoil surface is and the degree of Bezier curve is , then is a matrix, is a matrix, and is a matrix. Two columns of the matrix pertain to the - and -coordinates of the predetermined data on the airfoil surface. Equation (20) can be rewritten as If , the matrix will be a square matrix and it can be inverted. In such a case, (21) can be written as follows to find the matrix : However, the number of the airfoil surface data points is usually more than the number of control points. In such a case, there are more equations than unknowns and the matrix is no longer a square matrix. Hence it is required to convert it to a square matrix by multiplying both sides of (21) by the transpose of as follows: Thus, NACA 0015 and TsAGI “B” 12% airfoils produced by Bezier curve with and and their comparison with conventional NACA 0015 and TsAGI “B” 12% airfoils are shown in Figures 8 and 9, respectively. There is an excellent agreement between two airfoils in each figure.
The predetermined data for the NACA airfoils can be extracted from, for example, the software JavaFoil [37] which is based on the analytical NACA formulations.
NACA 00xx Symmetric Airfoils. Since the maximum thickness of a NACA 00xx symmetric airfoil will be considered as a design variable, the equation for generating such airfoils is given as follows: where is coordinates along the chord of the airfoil, from to ( is the chord length and is assumed equal to ), is the thickness coordinates above and below the line extending along the length of the airfoil, and is maximum thickness of the airfoil in percentage of chord (i.e., in a %15 thick airfoil would be 0.15). Equation (25) can be used to find the -coordinates of a NACA 00xx symmetric airfoil by knowing the values for and . As will be shown, the maximum thickness of such airfoils will also be considered as a design variable. By optimizing the thickness, the optimal shape for such airfoils will be obtained. This kind of optimization problem, however, is not comprehensive and produces the optimal NACA 00xx symmetric airfoils only. In summary, three kinds of design variable will be considered in this paper for airfoil shape optimization which are grid points on a given airfoil surface extracted from, say, the software JavaFoil, the Bezier curve control points, and the maximum thickness of NACA 00xx symmetric airfoils.
5. Shape Optimization
Different objective functions may be considered for the aerodynamics shape optimization including maximizing the lift-drag ratio, maximizing the lift, and minimizing the drag. In the framework of this paper, the shape optimization problem will be to infer the shape an airfoil should have so that the pressure distribution on the airfoil surface matches a prescribed one (an inverse problem). In inverse design problem, the desired pressure distribution of the target design may be specified a priori.
Design Variable (DV). Here the airfoil grid points, the Bezier curve control points, and the maximum thickness of NACA00xx airfoils are considered as design variables. Therefore, one has the following: Case 1: the airfoil grid points as design variable (see Figure 13); Case 2: the Bezier curve control points as design variable; Case 3: the maximum thickness of NACA00xx airfoils.
Case 1. The mathematical expression for the objective function considered for Case 1 can be stated as where is the pressure at grid points on the airfoil surface and is the desirable pressure at grid points on the airfoil surface (Figure 10). The aim is to minimize and to reach the desirable pressure distribution by changing the position of the grid points on the airfoil surface. Since the -coordinates of the grid points can be constant during the optimization process, only the -coordinates of the grid points are considered as design variables. Two end points of airfoil, namely, leading edge and trailing edge , are fixed. Thus they are not considered as design variables.
Case 2. The mathematical expression for the objective function considered for Case 2 can be stated as where is the number of the predetermined data on each of the upper and the lower surfaces of the airfoil, is the pressure at point of the airfoil surface generated by the Bezier curve, and is the desirable pressure at point . Why does ? data points for the upper surface, data points for the lower surface, and the leading and the trailing edges for two surfaces are considered fixed. The aim is to minimize and to reach the desirable pressure distribution by changing the -position of the control points on each of the upper and the lower surfaces of the airfoil (see Figure 11). and , which are concerned with the leading edge and the trailing edge, respectively, are considered fixed for both upper and lower surfaces. Therefore, for the shape optimization problem with a Bezier curve of order 11, we have 2 × (11 − 2) = 18 design variables. For the shape optimization problem with a Bezier curve of order 7, we have 2 × (7 − 2) = 10 design variables. The reason for considering these two kinds of the Bezier curve is twofold:(1)to show that the optimization problem will be more successful if we have less number of design variable;(2)to have a very accurate and flexible representation of the airfoil shapes, a degree of at least 10 should be used.
Case 3. The airfoil surface is generated by the analytical NACA formula (25) and the maximum thickness is considered as the design variable. To show the accuracy of the sensitivity scheme, the upper and lower airfoil surfaces are generated separately and hence the design variables will be two maximum thicknesses in the upper and lower airfoil surfaces. As shown in Figure 12, if the indices 1 and 2 denote the location of maximum thickness on the upper and lower airfoil surfaces, respectively, then the mathematical expression for the objective function considered for Case 3 is as follows:
6. Sensitivity Analysis
Suppose we wish to calculate the sensitivity of pressure of nodes on the airfoil surface (see Figure 10), , to the -position of the nodes on the airfoil surface, . The sensitivity analysis can be performed by introducing small perturbations to the -coordinate of each point on the airfoil surface, individually. The grid generation and flow problem may be solved for this perturbed shape to obtain the new values for the pressure . Using these values for the pressure, the dependency of the pressure to the perturbation of the -position of points of coordinates , , can be evaluated. The finite difference method may be used to formulate these sensitivities as follows: where may be, say, . The term is the perturbation in the -position of points of coordinates , . Since the sensitivity of each pressure to each -position of points of coordinates is required, the computation of the sensitivity coefficients using this method requires additional solutions of the flow problem. Therefore, this method is only suitable when the number of points on the airfoil surface is small. Thus for the airfoil shape optimization problem, which demand a fine grid to obtain accurate results, the perturbation method using the finite difference method will be of high computation cost. In this paper, we will expand the novel method used in evaluating the sensitivity matrix in the shape optimization of heat transfer problems. As will be shown, it requires only one solution of the flow problem (at each iteration) to compute all sensitivity coefficients.
With regard to (26), the following equation can be written in order to calculate the Jacobian matrix where and . The expression in the above relation is called the Jacobian coefficient. In this case, the sensitivity matrix can be expanded as
Since the physical domain is mapped onto the computational one, the chain rule may be used to correlate variables in the two domains. Therefore, As pointed out before, the -coordinate of the grid points are considered fixed and they are not included in the design variables. Thus (32) is written here to derive the required relations for the sensitivity coefficients. By interchanging and , and and , and solving the derived equations for and , we finally obtain where is the Jacobian of the transformation. Using the finite difference method to discretize the equations in the computational domain, we can write appropriate algebraic approximations for all partial derivatives involved in the above equation. Therefore, which are based on the central and the forward differences. Equations (35) through (38) are employed to calculate the sensitivity coefficients in (31).
Bezier Control Points as Design Variables. With regard to (27) and considering the control points of the Bezier curve as design variable, we can write Using the chain rule, we can write where are the -coordinate of the predetermined grid points to be passed by the Bezier curve and are the -coordinate of Bezier control points whose number is equal to 18 (9 for each of the upper and lower surfaces). The term can be computed by the expressions derived for Case 1 (see (31)). The size of the matrix formed by the arrays is . The term can be easily evaluated by taking derivative of (19) with respect to the control points (noting that ). The control points may be renumbered so that , , and , , . The indices and denote the upper and lower surfaces, respectively. The direction of numbering is from right to left for the upper surface and from left to right for the lower surface. The reason for this renumbering is the compatibility with the grid point data reading (most of the airfoil data are in this format) as well as the pressure reading to compute the objective function (see (27)). However, we should note that the Bezier curve evaluation is from left to right for both the upper and lower surfaces. The size of the matrix formed by the arrays is . Because the upper and lower surfaces are constructed separately, the variation of of the upper surface with respect to the change in position of the lower surface control points as well as the variation of of the lower surface with respect to the change in position of the upper surface control points is zero.
Maximum Thickness as Design Variables. In a similar derivation to Case 1, the sensitivity matrix for Case 3 will be
7. Optimization Method
In this paper, two powerful optimization methods, namely, the conjugate gradient method and the quasi-Newton method will be used. For the airfoil grid points as a design variable (Case 1) both optimization methods will be employed. However, for the Bezier curve control points as a design variable (Case 2), only the quasi-Newton method will be used. For Case 3 (the maximum thickness of NACA00xx airfoils as a design variable), only the conjugate gradient method will be employed.
Conjugate Gradient Method. The conjugate gradient algorithm to obtain the optimal shape for the airfoil is as follows.(1)Specify the physical domain, the boundary conditions, the problem conditions such as Mach number and the angle of attack, and the desired airfoil surface pressure distribution.(2)Generate the boundary fitted grids using the grid generation methods described earlier.(3)Solve the direct flow problem of finding the pressure values at any grid points of the physical domain and hence the airfoil surface.(4)Using (26), compute the objective function .(5)If the value of the objective function obtained in step (4) is less than the specified stopping criterion, the optimization is finished. Otherwise, go to step (6).(6)Compute the sensitivity matrix from (31).(7)Compute the gradient direction from (30).(8)Compute the conjugation coefficient from the following equation (the Polak-Ribiere formula): For , set .(9)Compute the direction of descent from the following: (10)Compute the search step size from the following: (11)Evaluate the new -coordinates of the airfoil surface grid nodes as follows: (12)Set the next iteration and return to step (2).The above algorithm is for the airfoil grid points as a design variable (Case 1) only. The algorithm for Case 3 can be expressed in a similar way.
Quasi-Newton Method. Quasi-Newton method is another powerful optimization method used in this paper. In quasi-Newton method, the Hessian matrix (which is composed of the second partial derivatives) is replaced by an approximation of it. The approximation uses only the first partial derivatives. The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method is a quasi-Newton method for solving unconstrained nonlinear optimization. In the BFGS method, the Hessian matrix approximation, , is updated iteratively. The steps of BFGS method can be summarized as follows.(1)Specify the physical domain, the boundary conditions, the problem conditions such as Mach number and the angle of attack, and the desired airfoil surface pressure distribution.(2)Generate the boundary fitted grids using the grid generation methods described earlier.(3)Solve the direct flow problem of finding the pressure values at any grid points of the physical domain and hence the airfoil surface.(4)Using (26), compute the objective function .(5)If value of the objective function obtained in step (4) is less than the specified stopping criterion, the optimization is finished. Otherwise, go to step (6).(6)Compute the sensitivity matrix from (31).(7)Compute the gradient direction from (30).(8)The initial Hessian matrix approximation, , is taken as the identity matrix, namely, .(9)Set and the iteration number as .(10)Compute the search step size (from (44)) in the direction and set (11)Repeat the steps (2) to (7) with these new values of for the grid points -coordinates to calculate .(12)Update the Hessian matrix approximation as where (13)Set the new iteration number as and go to step (9).
8. Results
In this section, the results obtained for the shape optimization of an airfoil in the incompressible, irrotational, and inviscid flow under given boundary conditions are presented. Three kinds of the design variable (the airfoil grid points, the Bezier curve control points, and the maximum thickness of NACA 00xx airfoils) as well as two optimization methods (CG and BFGS) are considered. In all test cases in this paper which employ the Bezier curve, the number of predetermined airfoil data, , is set equal to the Bezier curve order, .
Test Case 1. In this test case, the airfoil surface is parameterized by a Bezier curve of order 11. The total number of the design variable is 18, namely, 9 design variables for each of the upper and lower surfaces. At first, two parametric curves for two surfaces (upper and lower) are obtained using 11 grid points and then a fine grid is generated to obtain accurate results. The data for Test Case 1 is given in Table 1. The comparison of the initial and optimal airfoil shapes and some magnified parts of them are shown in Figures 15, 20, and 26. In this test case, BFGS optimization method is employed.
The convergence of the objective function is shown in Figures 17, 23, 29, 37, and 39. The initial and minimum values for the objective function are approximately 3517433 and 2877116, respectively, which shows %18.2 reduction in objective function (see Figures 16, 21, and 27). The minimum value for the objective function takes place in iteration 14. The optimization time spent on the 1st iteration (which is equivalent to one direct flow solution) is 11 minutes and 43 seconds and the total optimization time for 30 iterations is 11 minutes and 46 seconds which shows the proposed sensitivity analysis efficiency. 30 iterations take only 3 seconds. The reason for the difference between the 1st iteration and the following ones is that the solution after the 1st iteration is a very good initial guess for the 2nd iteration and the direct solution converges quickly. In other words, what is a bit time consuming for the 1st iteration is the grid generation and stream function loops not the pressure calculation, sensitivity analysis, and optimization stages. Moreover, a fine grid (300 × 305) and a tolerance of 10−8 are used in the iterative loops which increase the computation time. The code is programmed in Fortran 77 using a Fortran compiler (Force 2.0) and the computations are run on a PC with Intel Pentium Dual 1.73 and 1 G RAM. All the computations in the test cases in this paper are performed using the above mentioned compiler and PC. Therefore, there is no need to repeat it in the following test cases.
Test Case 2. Test Case 2 is similar to Test Case 1 but with different specifications (see Figure 18). The data for this test case is given in Table 2.
The explanation is similar to Test Case 1. Thus only the results will be given in Table 3.
Now an optimal shape design problem using a Bezier curve of order 7 is given. As it will be shown, it decreases the objective function value much bigger than when using a Bezier curve of order 11 as there is less number of design variables (10 design variables for a Bezier curve of order 7 versus 18 design variables for a Bezier curve of order 11). However, as it will be shown, using a Bezier curve of order 7 is not comprehensive for all airfoil shapes and is suitable to NACA 00XX or similar airfoils only. In other words, it is not able to produce all airfoil shapes precisely.
Test Case 3 (using a Bezier curve of order 7; see Figure 24 and Table 4). The results are given in Table 5.
Although there is a decrease of %59 in objective function, the %59 approach from the initial shape to desired one is not seen (see Figure 28). This indicates that the solution of the inverse problem is not unique. The reason for this can be found in the trailing edge configuration for the initial, optimal, and desired airfoil shapes (Figure 28).
Although the results of using the Bezier curve of order 7 is very promising, its drawback is that it is restricted to the simple and symmetric airfoil shapes such as NACA 00xx. For other airfoil shapes, there can be seen some oscillations around the trailing edge (see Figure 30).
Moreover, the Bezier curve of order 7 fails to represent the NACA 00xx airfoils accurately. In other words, the Bezier curve of order 11 is the appropriate option to produce very accurate airfoils (Figure 31).
Test Case 4. In this test case, the airfoil surface is parameterized by grid points obtained from the analytical NACA formula (e.g., software JavaFoil). In this case, the number of the design variables is equal to the number of grid points minus three (one for leading edge and two for trailing edge). Therefore, we have an aerodynamic shape optimization problem with a high number of the design variables as we should have a fine grid to obtain sufficiently accurate results. It is known that the optimization process becomes more challenging by increasing the number of design variables. Hence we have a difficult shape optimization problem in Test Case 4. The data used for Test Case 4 is given in Table 6. A very fine grid (400 × 425) is used for the initial airfoil shape (NACA 0012). The number of the design variables is which is 425 – 3 = 422. Therefore, a time consuming optimization problem is expected. However, by using the sensitivity analysis scheme proposed in this paper, the total time for the optimization problem in Test Case 4 using both CG and BFGS optimization methods is about 46 minutes for 8000 iterations. The computation time for the 1st iteration is about 25 minutes. The comparison of the computation times for the 1st iteration and 8000 ones indicates again the efficiency of the sensitivity analysis scheme. The summary of the results is presented in Table 7. The comparison of the initial, optimal, and desired airfoil shapes is given in Figure 32. As can be seen in the figure, the variation of the shape is minute. The convergence of the objective function to a local minimum when both the CG and BFGS optimization methods are used as well as a comparison of them is shown in Figures 33, 34, and 35, respectively. The plots reveal the better performance of the BFGS method in minimizing the objective function.
In Test Cases 5 and 6 the maximum thickness of the NACA 00xx is considered as a design variable. As mentioned previously, the conjugate gradient method is employed as the optimization method.
Test Case 5. The data for the problem including the conditions for the initial and desired airfoils are given in Table 8.
Figure 36 represents the comparison of the initial, optimal, and desired shapes for airfoils (also see Figures 22 and 38). The desired airfoil shape is a NACA0018 at conditions stated in Table 8. As can be seen, this shape is shown by a black color line. The optimization process is started by a NACA0011 as an initial shape which is shown by a red color line. The optimal shape (shown by a blue color line), which is obtained by the conjugate gradient method, is in an excellent agreement (full matching) with the desired one. The objective function variation is shown in Figure 37. The initial and final values for the objective function are about 773597.45 and 15.48, respectively, which reveals an approximately %100 reduction in the objective function within 31 iterations. The total time for the optimization (for 31 iterations) is 2 minutes and 14 seconds. The tolerance used in iterative steps in the program is 10−8. Although such a tolerance value increases the computation time, it enhances the accuracy of the results. If the -components of the maximum thickness in upper and lower airfoil surfaces are denoted by and , respectively, the value of the pressure for these two locations for initial, optimal, and desired shapes are reported in Tables 10 and 13. The difference values show the validity of the shape optimization process.
Test Case 6. The data for the problem is given in Table 11. The explanation for the results is similar to Test Case 5 (see Tables 9 and 12).
9. Adjoint Method
As pointed out previously, for the aerodynamic shape optimization problems requiring a large number of design variables, the use of finite difference method to evaluate the gradient by introducing a small perturbation to each design variable separately and then solving the flow problem is of very high computational cost, because it requires a number of additional flow solutions equal to the number of design variables. For optimal shape design problems with a high number of design variables, the adjoint method [4] can compute the gradients of objective function much faster than the finite difference method.
The aerodynamic shape optimization problem of interest here can be expressed asThe objective function and the governing equation depend on the flow variables and the geometry design variable : The derivative of the objective function with respect to the design variables can be expressed as which states that a change in the objective function is due to a combination of a variation in the flow solution and a variation in the design variable (change in geometry) . In a similar way, we have If the sensitivity analysis is performed using (51) and (52), the problem is referred to as the “primal problem.” Solving the primal problem comes with the same difficulties as we encounter with use of the finite difference method. It requires the additional flow solutions proportional to the number of the design variables . Therefore, the adjoint method comes to the picture by introducing a vector of Lagrange multipliers . By adding (52) as a constraint to the sensitivity equation (51), we obtain Rearranging the terms inside (53), we get If then (54) reduces to Equation (55) is the adjoint equation and the vector is the adjoint variables. Equations (55) and (56) are referred to as the “dual problem.” The adjoint equation is a linear system and can be solved to obtain . Then the determined can be substituted into (56) to obtain the gradient of the objective function. It can be seen that the gradient of the objective function can be determined without the need for additional flow solutions. The computational cost of solving the adjoint equation is comparable to that of solving the flow equation. Therefore, the computational cost of evaluating the objective function gradient is roughly equal to the computational cost of two flow equation solutions, independent of the number of design variables [38–41].
From the accuracy of the derivatives view point, the finite difference method (based on the perturbation scheme) is compared to the adjoint method [42–44]. The comparison shows a very good agreement between two methods. Therefore, our aim here is to compare our novel shape sensitivity method to the adjoint method from the efficiency view point only. As mentioned above, the computational cost of solving the adjoint equation is comparable to that of solving the flow equation whereas the computational cost of our novel method is comparable to that of computation of an algebraic expression for arrays of a matrix. As seen in Test Case 4, the computation time for iterations 2 to 8000 is 46 – 25 = 21 minutes (about 7 iterations per second) which reveals the efficiency of the proposed sensitivity analysis.
10. Conclusion
This paper addressed the aerodynamic shape optimization for an airfoil in an irrotational and incompressible flow governed by Laplace equation using a type of the elliptic grid generation (O-type), a novel and very efficient sensitivity analysis method, and the conjugate gradient and BFGS optimization methods. The airfoil was parameterized using the grid points and the Bezier curve. Three different types of design variable were considered: the grid points, the Bezier curve control points, and the maximum thickness of NACA00xx airfoils. It was represented that the use of the Bezier curve significantly improves the optimization performance to reach the optimal shape. The results obtained in test cases presented in this paper show that the proposed sensitivity analysis method reduces the computation cost even for large number of the design variables (Test Case 4) and confirm accuracy and efficiency of the proposed shape optimization algorithm.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.