Next Article in Journal
The Clash between CLIL and TELL: Effects and Potential Solutions of Adapting TELL for Online CLIL Teaching
Next Article in Special Issue
Design and Control of a Linear Rotary Electro-Hydraulic Servo Drive Unit
Previous Article in Journal
Feature Analysis of Predictors Affecting the Nidus Obliteration of Linear Accelerator-Based Radiosurgery for Arteriovenous Malformations Using Explainable Predictive Modeling
Previous Article in Special Issue
A Multi-Objective Optimization Method of a Mobile Robot Milling System Construction for Large Cabins
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Robust Algorithm Software for NACA 4-Digit Airfoil Shape Optimization Using the Adjoint Method

by
Naser Tanabi
,
Agesinaldo Matos Silva, Jr.
,
Marcosiris Amorim Oliveira Pessoa
and
Marcos Sales Guerra Tsuzuki
*
Computational Geometry Laboratory, Departamento de Engenharia Mecatrônica e de Sistemas Mecânicos, Escola Politécnica da Universidade de São Paulo, São Paulo CEP 05508-030, Brazil
*
Author to whom correspondence should be addressed.
Appl. Sci. 2023, 13(7), 4269; https://doi.org/10.3390/app13074269
Submission received: 27 December 2022 / Revised: 24 March 2023 / Accepted: 26 March 2023 / Published: 28 March 2023

Abstract

:
Optimizing the aerodynamic shape of an airfoil is a critical concern in the aviation industry. The introduction of flexible airfoils has allowed the shape of the airfoil to vary, depending on the flight conditions. Therefore, in this study, we propose an algorithm that is capable of robustly optimizing the shape of the airfoil based on variable parameters of the airfoil and flight conditions. The proposed algorithm can be understood as an optimization method, which employs the adjoint method, a powerful tool for estimating the sensitivity of the model output to the input in numerous studies. From an aerodynamic perspective, the development of shape geometry is a crucial step in airfoil development. The study used NACA-4 digit airfoils as input for the initial assumption and the range of shape change. The optimal shape was found using the proposed algorithm by defining one NACA profile as the initial value and another NACA profile as the limit for the optimized shape, considering the aerodynamic coefficients and flight conditions. However, morphing airfoils have certain deformation limitations. As an innovation in the algorithm, bounds were defined for the shape change during optimization so that the result can be constructed within the capabilities of the morphing wing. These bounds can be adjusted (depending on the capabilities of the airfoils). To validate the proposed algorithm, the study compared it with a previous flow solver for the same airfoil.

1. Introduction

Determining the geometric shape of an aircraft’s wings is crucial for optimizing its aerodynamics and performance. Adjusting geometric constraints can improve various aspects, such as aerodynamic performance, flight distance, and flight time. The first steps in wing design involve defining the wing’s geometry and selecting an appropriate airfoil type based on specific parameters. Special mission requirements may require designing and manufacturing the necessary parts and equipment to meet different operating conditions. The wing and airfoil are critical components, and numerical optimization techniques can be divided into two main categories: gradient-free and gradient-based methods. When dealing with large-scale aerodynamic optimization problems with multiple variables, the adjoint method is the most effective approach, particularly when the number of design variables exceeds 100 [1,2].
The focus of this paper is on robust aerodynamic shape optimization; we present an algorithm that differs from other similar articles. While other articles rely on commercial software or algorithms that do not consider flexible airfoil capabilities, this study takes into account the airfoil’s global geometry instead of solely modifying the mesh points of its boundary without regard to its flexibility. Building upon previous research by Tanabi et al. [3], who utilized the pressure distribution around the airfoil to determine its shape [3], this work uses design parameters and flight conditions to determine the airfoil’s shape, which is a more adaptable and versatile approach.
The primary objective of this study is to develop an algorithm that considers critical parameters, such as chord length, maximum thickness, and maximum camber, for flexible NACA-4-digit airfoils. The algorithm employs a unique reverse design approach to determine the shape of the NACA-4 digit airfoil based on the surrounding pressure distribution. During each iteration of the algorithm, it checks if the pressure distribution can generate NACA-4 digits while imposing constraints to ensure that the variables remain within acceptable ranges.
In this study, we conducted a comprehensive review of relevant literature and clarified the adjoint equation and optimization algorithm. The flow solver was validated by comparing its results to those obtained through numerical solutions. Finally, we present the results of the multipurpose optimization approach applied to an airfoil in a given regime, followed by a detailed discussion of the proposed algorithm’s performance.

2. Bibliographic Review

The introduction of the gradient-based CFD method by Jameson [4] in 1988 marked the first application of adjoint equations to transient flow [4]. They are also among the most prominent developers of this method and have utilized adjoint equations in various computational fluid dynamics test cases, including the optimization of aerodynamic shapes for Eulerian flows that govern the Navier–Stokes equations. Lions was the first person to use optimization to solve differential equations [5], while adjoint equations were originally developed for shape optimization in fluid dynamics to solve a defined control problem [5,6]. The coupling of perturbation equations with this method was achieved by integrating the Reynolds averaged Navier–Stokes (RANS) equations with the Spalart–Allmaras single equation turbulence model [7,8]. Several researchers, including Giles et al. [9], utilized the adjoint optimization method in fluid dynamics and Eulerian equations and conducted significant research in this field [9,10]. Later, Giles and Pierce [11] used Green’s functions for the adjoint equations [11]. To design an airfoil, they used the adjoint method under incompressible non-viscous conditions in the far-field region [12], while Xie [13] used the adjoint equations by considering the angle of attack based on the boundary conditions of the distant region [13].
In the multi-objective optimization field, various types of research studies have been conducted. Braibant and Fleury [14] presented a method for shape optimization utilizing B-splines, which are piecewise polynomial curves that can represent intricate shapes, making them valuable in engineering designs [14]. Poles et al. [15] investigated the effect of initial population sampling on the convergence of multi-objective genetic algorithms (MOGAs) [15]. The authors conducted experiments using different techniques for initial population sampling and evaluated their impact on the performance of MOGA. Additionally, they introduced a new evolutionary algorithm for multi-objective optimization, MOGA-II, and tested its performance on single-objective optimization problems, comparing it to two successful algorithms: differential evolution and a standard evolutionary algorithm. Mariotti et al. [16] conducted experiments using particle image velocimetry (PIV) and discovered that incorporating contoured cavities on the upper surface of the diffuser helps to reduce the extent of the flow separation region, thus enhancing the diffuser’s efficiency. They also presented a study that aimed to improve diffuser efficiency using multiple local recirculation regions (MLRR) [17].
The Aerodynamic Shape Optimization Design Group is currently researching optimization methods and developing new techniques and algorithms. In recent studies, they have employed an adjoint-based solver to compute sensitivity derivatives. This approach is demonstrated in two articles: Telidetzki et al. [18] used the B-Spline method to adjust volume parameters, and calculated sensitivity derivatives using an adjoint equation solver and a jet stream solver. This approach was used to obtain an optimized shape for the NACA 0012 airfoil [18].
In a study by Amoignon et al. [19], three methods for the shape optimization of airfoils were compared [19]. The first two methods used free-form deformation, and the third method used the radial basis function (RBF) coupled with the adjoint method. Sensitivity derivatives were obtained using an adjoint-based solver with an unstable flow solver and a finite difference method. The drag coefficients of the NACA 0012 and RAE 2822 airfoils were reduced under free deformation conditions. In an article by Carrier et al. [20], the Bezier curve technique was used to modify the shapes of the NACA 0012 and RAE 2822 airfoils for shape optimization [20]. The sensitivity derivatives were calculated using a flow solver and an adjoint solver. The conjugate gradient method was used as an optimization algorithm in this study. Carrier et al. [20] achieved a considerable reduction in the pressure drop ratios of NACA 0012 and RAE 2822 using this method [20].
Adjoint optimization has gained popularity, especially in modern industries. In a study by Othmer et al. [21], the method was applied to analyze a transient flow through 2D and 3D ducts using a solver that combined a finite-volume method with continuous-adjoint optimization. Another project was carried out in collaboration with Volkswagen Automotive [22], where the group utilized adjoint equations to improve the aerodynamic shape of various components such as cockpit air ducts, engine components, propulsion, exhaust flow, and the car body. One particular area of interest was the flow around the side mirrors, which involved modifying the mirror surfaces. In 2004, the adjoint method was also used to investigate a multipurpose aerodynamic shape optimization method with a slotted flap [23]. The study utilized a single and multipurpose algorithm and the Newton–Krilov gradient-based method, with twelve moving points on the surface of an airfoil and eight points on the surface of the flap.
In a study by Schramm et al. [24], an adjoint optimization algorithm was employed using the finite difference method to optimize the shape of the front edge (leading edge flaps) of a slotted airfoil, resulting in improved aerodynamic coefficients. Furthermore, the optimized shape was experimentally validated in a wind tunnel. In a separate study, Rashad and Zingg [25] utilized the discrete adjoint method to optimize the aerodynamic shape of an airfoil by simulating natural laminar flow, which included a transient turbulent current under the RANS method. The adjoint method has become a prevalent tool for shape optimization; Elham and van Tooren [26] used this technique in OpenFoam software to optimize the shape of an airfoil [26].

3. Methodology

This section presents the proposed algorithm, which consists of several components. First, the NACA four-digit airfoil is introduced. Then, the basic concepts of the adjoint equation are explained, taking into account the flow field, the boundary geometry, and the cost function. The algorithm developed can handle both viscous and non-viscous flows, and in this section, the adjoint equation for viscous flow is described. Finally, the optimization algorithm is presented.

3.1. Geometry Generation of the NACA Four-Digit Airfoil

The NACA airfoil was initially developed for aircraft, but it has also found applications in wind turbines. In the literature, the four- and five-digit types of NACA (National Advisory Committee for Aeronautics) have been extensively investigated. The coordinates of the airfoil ( x , y ) in a plane ( x y ) are not dimensionless with respect to the chord. Airfoils with NACA four-digit designations have piecewise parabolic mean camber lines y c and thickness distributions y t that are consistent across designs. The thickness distribution is given by
y t = ± 5 t c 0.2969 x 0.1260 x 0.3516 x 2 + 0.2843 x 3 0.1015 x 4
where t is the thickness ratio and c is the chord length of the airfoil. The mean camber line y c comprises two sections, which are defined by different parabolic curves, as in
y c = m p 2 ( 2 p x x 2 ) 0 x < p m ( 1 p ) 2 ( ( 1 2 p ) + 2 p x x 2 ) p x 1 ,
where p is the maximum camber position along the chord length and m is the camber ratio. In other words, the complete set of the coordinates ( x , y ) of the NACA airfoil can be described by two curves, i.e., the upper and lower surfaces, which are defined by
y upper = y c + y t cos ( θ ) x upper = x y t sin ( θ ) y lower = y c y t cos ( θ ) x lower = x + y t sin ( θ ) ,
where,
θ = t a n 1 d y c d x .
To summarize, the first digit of a NACA airfoil indicates the maximum mean camber, the second digit indicates the position of the maximum camber, and the last two digits indicate the maximum thickness. The NACA4418 airfoil selected as an initial condition for the optimization process described in the following is an asymmetric airfoil with a maximum mean camber of 4% located at a distance of 0.4 c and a maximum thickness of 18%, which is t = 0.18 .

3.2. Adjoint Equation

The progress made in design procedures relies on the cost function (I). Parameters such as the pressure distribution on an object or the lift-to-drag ratio in reverse problems may be used to evaluate the cost function. When dealing with the flow around aerodynamic structures or wings, the cost function (I) represents aerodynamic features and involves the flow field w and the boundary geometry F . Therefore, the cost function may be viewed as a function of these two variables.
I = I ( w , F ) .
As a result, modifying the boundary F can cause a corresponding alteration in the cost function
δ I = I w T δ w + I F T δ F .
The initial term pertains to alterations in the cost function with the flow variables, whereas the subsequent term pertains to modifications in boundary geometry variables. When using control theory, flow field equations should not generate multiple solutions, ensuring the removal of δ w from Equation (6). If we denote the flow equations as R, their general forms can be categorized into two types, based on the flow and geometric variables
R = R ( w , F ) .
As a result, modifying the boundary F leads to
δ R = R w δ w + R F δ F = 0 .
When applying the finite difference method, one can obtain δ w from Equation (8) and insert it into Equation (6). However, this approach is not practical in complex and multivariate problems, such as fluid flow fields, and can lead to issues with convergence and obtaining the correct solution. The adjoint method offers an alternative solution by removing δ w -dependent terms from Equation (6). This can be achieved by defining the Lagrangian coefficients ψ and assuming that δ R = 0 . Then, δ I can be expressed as
δ I = δ I ψ T δ R δ I = I w T δ w + I F T δ F ψ T R w δ w + R F δ F .
The above sentences can be rearranged as
δ I = I w T ψ T R w δ w + I F T ψ T R F δ F .
The arbitrariness of the Lagrangian coefficients ψ allows for the first component of Equation (10) to be zeroed, as
I w T ψ T R w = 0 .
By eliminating the dependencies of the changes of the cost function δ I on the flow variable δ w , only the geometric boundary variable remains. Equation (11) is referred to as the adjoint equation, where the Lagrangian coefficients ψ are the unknowns. When this equation is applied to Equation (10), G , the gradient of the cost function (with respect to the boundary geometric variables) is obtained as follows
G = δ I δ F = I F T ψ T R F .

3.3. Adjoint Equation for Viscous Flow

The equation that describes the overall flow in a domain Ω can be stated as
w t + f i x i = f v i x i in Ω .
The vector w corresponds to the state variables, and f i and f v i are the viscous flux and non-viscous flows for a given dimension x i , respectively. In the case of a two-dimensional flow field that is defined by x 1 = x , x 2 = y , f 1 = f , and f 2 = g , the resulting expression is obtained as
w t + f x + g y = f v x + g v y in Ω ,
where the vectors g and g v correspond to the viscous and non-viscous fluxes in the y directions, and they are defined as follows
w = ρ ρ u ρ v ρ E , f = ρ u ρ u 2 + p ρ u v ρ u H , g = ρ v ρ v u ρ v 2 + p ρ v H , f v = 0 σ x x σ y x u σ x x + v σ x y + k T x , g v = 0 σ x y σ y y u σ x y + v σ y y + k T y .
The steady-state non-viscous flow equation by neglecting the viscous sentences of Equation (14) can be written as
R = f x + g y = 0 in Ω .
The relationship can be defined as R and can be substituted for the adjoint Equation (11), resulting in
I w T ψ T R w = 0 in Ω , I w T ψ T w f x + g y = 0 in Ω , I w T ψ T x f w + y g w = 0 in Ω .
Integrating the volume will yield
Ω I w T Ω ψ T x f w + y g w d Ω = 0 in Ω .
The second integral, using the partial integrals, is as follows
Ω ψ T x f w + y g w d Ω = B ψ T f w n x + g w n y d B Ω ψ T x f w + ψ T y g w d Ω .
The first volume integral of Equation (19) will also be commutable in the form of a cost function. For example, it assumes that the cost function for achieving the optimal pressure on the solid boundary (inverse problem) is defined as
I = 1 2 B p d i s p d i s * 2 d B .
In this case, since this cost function is defined on boundary B, its volume integral will be zero, as
Ω I w T = 0 in Ω .
After substituting Equations (21) and (19) into the integral adjoint Equation (18), the resulting expression is obtained as
B ψ T f w n x + g w n y d B Ω ψ T x f w + ψ T y g w d Ω = 0 .
The condition for this equation to be zero is that the integers of each integral are zero, separately, as
ψ T x f w + ψ T y g w = 0 in Ω ,
ψ T f w n x + g w n y = 0 in B .
The adjoint Equation (23) is formulated for the entire field with boundary conditions that satisfy Equation (24). To solve the adjoint Equation (23), an artificial term can be added to transform the equation into a time-dependent equation, as
ψ T t + ψ T x f w + ψ T y g w = 0 in Ω .

3.4. Optimization Algorithm

Based on the discussion of the adjoint equation formation for a fluid problem, the optimization algorithm can be conceptualized in the following manner:
  • Step 1: The initial step involves defining the initial parameters for the algorithm, which initiates the primary optimization loop.
  • Step 2: The flow field equations are solved to obtain a reliable solution. AUSM+ [27] was employed to compute the flow field variables that will be utilized in subsequent steps.
  • Step 3: The adjoint equation is formulated and solved by the algorithm. The algorithm reconstructs the adjoint equation using the flow field variables.
  • Step 4: The boundary geometry is deformed in the direction of the maximum steepest descent, adhering to the NACA standard. After observing the trend of movement of the boundary geometry, the shape of the airfoil is altered to approach the optimal condition.
  • The field mesh is reproduced, and the algorithm returns to Step 1 to attain the minimum-cost function. In each iteration, the shape undergoes modifications, and a new mesh is produced at the field boundary.
The loop of the algorithm aimed at achieving the optimized shape is depicted in Figure 1. To identify the path of deformation, the relation given by Equation (12) can be employed to enhance the form, as
δ F = λ G ,
where λ is the step size. The following equation provides the search direction toward the optimal value using the steepest descent method
S = G
Determining the step size is also necessary to move in this direction
λ = e b ln 10 e a ln 10 = e ( b a ) ln 10 ,
where a and b are the results of
a = 2 if S = 0 log | S | if S 0 , b = 2 if F = 0 log | F | if F 0 .
To maximize the ratio of the lift-to-drag coefficient, one can combine the objective functions of maximizing the lift coefficient and minimizing the drag coefficient. For this purpose, the coefficients of the general objective function must be calculated as
I = 1 C l C l max 2 if C d C d min 1 2 1 C l C l max 2 + 1 2 1 C d C d min 2 if C d > C d min ,
where C d min is the target value of the drag coefficient and C l max is the target value of the lift coefficient. As a result, the objective function changes, as
δ I = 2 1 C l C l max δ C l C l max if C d C d min 1 C l C l max δ C l C l max 1 C d C d m i n δ C d C d m i n if C d > C d m i n .
In the optimization process, limitations are set by taking into account the varying abilities of morphing airfoils compared to other airfoils. The range of changes allowed in the determining components of the flap geometry during the optimization path are presented in Table 1.

4. Results

The study case was chosen to assess the accuracy of the proposed algorithm in transient and high-speed flows. The optimization was performed under viscous flow conditions. The adjoint sensibility for 10 different grid generations in one optimization loop, the parameter change rate limitation, the grid convergence analysis for three grid generations by the Roache method, and the results of the multipurpose optimization algorithm are presented.

4.1. Study Case Selection

In a previous study, Tanabi et al. [3] evaluated an algorithm under laminar flow conditions [3]. The authors of the current study selected a different sample to evaluate the proposed algorithm under non-viscous conditions. They investigated the performance in transient flow by simulating the flow around a NACA 0012 airfoil with free-flow conditions of M = 0.85 and α = 1 . These parameters were chosen to assess the accuracy of the proposed algorithm in transient and high-speed flows. The simulation used a circular grid with dimensions of 100 × 70 .
A circular grid with dimensions of 100 × 70 was used to simulate the flow around the airfoil. Figure 2a,b show the lines of the pressure coefficient and Mach number, respectively. The simulation accurately predicted the two shock waves on the surface and below the surface, with an even lower surface shock wave captured, possibly due to the finer size of the grid compared to that used by Liou and Steffen Jr. [28]. The study compared the results of two solver methods, Roe and AUSM, for the pressure coefficient on the airfoil surface shown in Figure 2c, with the numerical results obtained in previous research.

4.2. Adjoint Sensitivity and Parameter Changing Range

To validate the adjoint sensitivity, the same input used in the example presented in the paper was applied to different grid conditions. The behavior of the adjoint sensitivity was observed in 10 different runs, which produced similar results as shown in Figure 3. However, additional discussion on this topic may divert the attention from the primary focus of the paper and overwhelm one with additional information. Therefore, we did not address this topic in the paper.
In the optimization process, constraints are established based on the morphing capabilities of the airfoils, which may vary from other types of airfoils. The allowable ranges for modifications to the determining components of the flap geometry in the optimization process are outlined in Table 2.

4.3. Grid Convergence Analysis

The order of convergence can be determined using the method proposed by Roache [29], based on the results obtained from testing the developed algorithm with different mesh generations. Figure 4 illustrates the residual of the adjoint variable for the main loop with 2400 iterations using 3 different mesh numbers of 50 × 35 , 100 × 70 , and 200 × 140 . The method provides the means to calculate the order of convergence using
p G C I = ln f 3 f 2 f 2 f 1 / ln ( r ) ,
where r is the constant grid refinement ratio, and f 1 , f 2 , and f 3 are the lowest residuals of the adjoint variables from cases 1, 2, and 3, respectively. For these test cases, the grid refinement ratio was r = 2 and the resulting lowest logarithmic residuals were log ( f 1 ) = 3.931 , log ( f 2 ) = 3.810 , and log ( f 3 ) = 3.796 . The order of convergence of p G C I = 2.848 was obtained for the analyzed data. The method provides the grid convergence index for a given relative difference between two cases using the following expression
G C I = F S · | ε | r p G C I 1 × 100
where ε is the relative error between two test cases as in ε i , j = f j f i / f i and F S is the safety factor. Since three grids were used to estimate p G C I , a safety factor of F S = 1.25 was used. The resulting CGI between Case 1 and 2 was G C I 1 , 2 = 0.6566 and between Case 2 and 3 was G C I 2 , 3 = 4.8875 .
To confirm the accuracy of the solution, the following relationship can be used with the three grids
G C I 3 , 2 r p × G C I 2 , 1 1 .
In our case, it is 1.0337 1 .

4.4. Multipurpose Optimization of Airfoil

The optimization process involves selecting the NACA4418 airfoil as the primary airfoil and using it to determine the final shape of the airfoil by minimizing the total drag coefficient and maximizing the lift coefficient based on design parameters and the primary airfoil’s shape. Optimization is performed under viscous flow conditions as specified in Table 3. The profile design variables include the maximum thickness, maximum camber, and maximum camber position, and their values are adjusted to achieve the intended purpose and obtain the final shape of the profile. The lift and drag values are deliberately chosen in such a way that they cannot be achieved under the defined flight conditions. The purpose of this work is to show the ability of the algorithm in defining the design limits of the shape of the airfoil. which moves toward optimizing the shape as much as possible; after reaching the defined design limits, it continues to optimize so that it can display the most ideal and closest possible shape to the output, according to the optimization goals.
Table 4 indicates that although the airfoil thickness and maximum camber decreased, and the camber location changed toward the end of the airfoil, only the drag coefficient was optimized. Therefore, it is not possible to design an airfoil that optimizes both conditions to meet them simultaneously under the specified flow conditions. In Figure 5, the convergence of the cost function, lift-to-drag ratio, and aerodynamic coefficients is illustrated. The graphs show that all parameters have converged to an acceptable value after 20 iterations. The aerodynamic coefficients shown in the figures include the viscous drag coefficient C d v , the pressure drag coefficient C d p , the lift coefficient C l , the total drag coefficient C d , and the lift-to-drag ratio C l / C d . Although both optimization goals were not achieved, the changes in all design parameters and aerodynamic coefficients have become insignificant.
The authors deliberately chose input values for maximum lift and minimum drag coefficients that were unattainable due to design limitations to evaluate the algorithm’s performance. The results show that the algorithm increases the lift-to-drag ratio by reaching the maximum lift target value while minimizing drag, but both goals cannot be achieved simultaneously due to design constraints. Figure 6 illustrates that after fifteen cycles, the maximum thickness t, the bend location m, and the maximum bend p exhibit insignificant changes.
The modifications in the pressure coefficient C p around the airfoil generated by the multipurpose optimization algorithm are shown in Figure 7. The lift coefficient is represented by the enclosed area in this shape. When comparing the distribution of the initial and final airfoil pressure coefficients surrounding the airfoil, it can be seen that the lift coefficient has increased. The shape of the airfoil deformation in the pattern to reach the final shape is illustrated in Figure 8. The algorithm attenuates the thickness until near the design limit (to reach the maximum lift coefficient). After that, it changes other design parameters, such as the maximum camber and position, to achieve the defined goal. The Mach contour around the airfoil is presented in Figure 9. There was a shock in the middle of the airfoil in the initial shape; after optimizing the shape, the shock appeared. The streamlines around the airfoil in the initial and final conditions are shown in Figure 10. The eddies were reduced after optimization, and the separation rate decreased in the airfoil tail, which is evident in the streamlines in the airfoil tail.

5. Conclusions

This paper presents an airfoil shape optimization method based on the adjoint approach. The non-viscous flow computations were carried out using Euler equations, while Navier–Stokes equations were used for viscous flow calculations. The AUSM+ scheme was utilized to calculate the flow coefficients. The adjoint method was used to apply the optimization equations to the flow equations, resulting in an algorithm for airfoil shape improvement.
As an innovation, multi-objective optimization was performed to evaluate the software performance, in which the lift and drag coefficients were determined as optimization goals. In the issue mentioned above, viscous flow was used. The airfoil design parameters were considered as the maximum thickness, maximum camber, and maximum camber position.
In the multi-objective optimization problem, it was observed that the total drag coefficient achieved the set value; however, the final lift coefficient had a significant difference from the set value. The outcome suggests that both sets of target values could not be achieved simultaneously under the flow conditions mentioned in the problem. Hence, it is not possible to find every coefficient value under a single flight condition. Consequently, to produce specific lift or drag coefficients, one must modify both the airfoil shape and flow conditions.

Author Contributions

Conceptualization, N.T.; methodology, N.T.; software, N.T.; validation, N.T.; formal analysis, N.T.; investigation, N.T.; resources, M.S.G.T.; data curation, N.T.; writing—original draft preparation, N.T. and A.M.S.J.; writing—review and editing, N.T., A.M.S.J. and M.S.G.T.; visualization, N.T.; supervision, M.S.G.T.; project administration, M.S.G.T. and M.A.O.P.; funding acquisition, M.S.G.T. All authors have read and agreed to the published version of the manuscript.

Funding

N. Tanabi was supported by Petrobras/ANP/FUSP. M. S. G. Tsuzuki is partially supported by CNPq (Grant 305.959/2016–6). The paper has the support by CAPES/PROAP-Grant 817.757/38.860.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

α Angle of Attack
λ Step Size of Steepest Descent
B Airfoil’s Boundary Geometry
F Boundary Geometric Variables of Cost Function
G Gradient of the Cost Function
S Search Direction of Steepest Descent Method
Ω Open Domain
ψ Lagrangian Coefficient
ε Relative Error
aDirection-based Step Size Coefficient
bGeometry-based Step Size Coefficient
cChord Length
C d p Pressure Drag Coefficient
C d v Viscous Drag Coefficient
C d Total Drag Coefficient
C l Lift Coefficient
C d min Target Drag Coefficient
C l max Target Lift Coefficient
fViscous Flux
F S Factor of Safety
f v Non-viscous Flux
G C I Grid Convergence Index
ICost Function
MMach Number
mMaximum Camber
n x X-component of Normal to B
n y Y-component of Normal to B
pMaximum Camber Position
p d i s * Target Pressure Distribution on Airfoil
p d i s Pressure Distribution on Airfoil
p G C I Order of Convergence
RFlow Equations
rGrid Refinement Ratio
R e       Reynolds Number
tMaximum Thickness
wFlow variables of Cost Function
y c Mean Camber Line
y t Thickness distribution

References

  1. Pulliam, T.; Nemec, M.; Holst, T.; Zingg, D. Comparison of evolutionary (genetic) algorithm and adjoint methods for multi-objective viscous airfoil optimizations. In Proceedings of the 41st Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 6–9 January 2003; p. 298. [Google Scholar]
  2. Lyu, Z.; Kenway, G.K.W.; Martins, J.R.R.A. Aerodynamic Shape Optimization Investigations of the Common Research Model Wing Benchmark. AIAA J. 2015, 53, 968–985. [Google Scholar] [CrossRef] [Green Version]
  3. Tanabi, N.; Barari, A.; Vakilipour, S.; Tsuzuki, M.S.G. Inverse Designing Airfoil Aerodynamics in Compressible Flow by Target Pressure Distribution. In Proceedings of the ICGG 2020-Proceedings of the 19th International Conference on Geometry and Graphics, São Paulo, Brazil, 9–13 August 2021; Springer: Berlin/Heidelberg, Germany, 2021; pp. 308–319. [Google Scholar]
  4. Jameson, A. Aerodynamic design via control theory. J. Sci. Comput. 1988, 3, 233–260. [Google Scholar] [CrossRef] [Green Version]
  5. Lions, J.L. Optimal Control of Systems Governed by Partial Differential Equations; Springer: Berlin/Heidelberg, Germany, 1971; Volume 170. [Google Scholar]
  6. Mohammadi, B.; Pironneau, O. Applied Shape Optimization for Fluids; Oxford University Press: Oxford, UK, 2010. [Google Scholar]
  7. Orovio-Kiachagias, E.M.; Asouti, V.G.; Giannakoglou, K.C.; Gkagkas, K.; Shimokawa, S.; Itakura, E. Multi-point aerodynamic shape optimization of cars based on continuous adjoint. Struct. Multidiscip. Optim. 2019, 59, 675–694. [Google Scholar] [CrossRef]
  8. Bueno-Orovio, A.; Castro, C.; Palacios, F.; Zuazua, E. Continuous adjoint approach for the Spalart-Allmaras model in aerodynamic optimization. AIAA J. 2012, 50, 631–646. [Google Scholar] [CrossRef]
  9. Giles, M.B.; Pierce, N.A.; Giles, M.; Pierce, N. Adjoint equations in CFD-Duality, boundary conditions and solution behaviour. In Proceedings of the 13th Computational Fluid Dynamics Conference, Snowmass Village, CO, USA, 29 June–2 July 1997; p. 1850. [Google Scholar]
  10. Giles, M.B.; Pierce, N.A. On the properties of solutions of the adjoint Euler equations. Numer. Methods Fluid Dyn. VI 1998, 1–16. [Google Scholar]
  11. Giles, M.B.; Pierce, N.A. Analytic adjoint solutions for the quasi-one-dimensional Euler equations. J. Fluid Mech. 2001, 426, 327–345. [Google Scholar] [CrossRef] [Green Version]
  12. Qiao, Z.; Yang, X.D.; Zhu, B. Numerical optimization design of wings by solving adjoint equations. Coordinates 2002, 3, 1. [Google Scholar]
  13. Xie, L. Gradient-Based Optimum Aerodynamic Design Using Adjoint Methods; Virginia Polytechnic Institute and State University: Blacksburg, VA, USA, 2002. [Google Scholar]
  14. Braibant, V.; Fleury, C. Shape Optimal-Design Using B-Splines. Comput. Methods Appl. Mech. Eng. 1984, 44, 247–267. [Google Scholar] [CrossRef]
  15. Poles, S.; Fu, Y.; Rigoni, E. The Effect of Initial Population Sampling on the Convergence of Multi-Objective Genetic Algorithms. In Proceedings of the Multiobjective Programming and Goal Programming: Theoretical Results and Practical Applications; Lecture Notes in Economics and Mathematical Systems; Barichard, V., Ehrgott, M., Gandibleux, X., Kindt, V.T., Eds.; Springer: Berlin/Heidelberg, Germany, 2009; Volume 618, p. 123. [Google Scholar]
  16. Mariotti, A.; Grozescu, A.N.; Buresti, G.; Salvetti, M.V. Separation control and efficiency improvement in a 2D diffuser by means of contoured cavities. Eur. J. Mech. B -Fluids 2013, 41, 138–149. [Google Scholar] [CrossRef]
  17. Mariotti, A.; Buresti, G.; Salvetti, M.V. Use of multiple local recirculations to increase the efficiency in diffusers. Eur. J. Mech. B-Fluids 2015, 50, 27–37. [Google Scholar] [CrossRef]
  18. Telidetzki, K.; Osusky, L.; Zingg, D.W. Application of jetstream to a suite of aerodynamic shape optimization problems. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA, 13–17 January 2014; p. 0571. [Google Scholar]
  19. Amoignon, O.; Hradil, J.; Navratil, J. Study of parameterizations in the project CEDESA. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA, 13–17 January 2014; p. 0570. [Google Scholar]
  20. Carrier, G.; Destarac, D.; Dumont, A.; Meheut, M.; Salah El Din, I.; Peter, J.; Ben Khelil, S.; Brezillon, J.; Pestana, M. Gradient-based aerodynamic optimization with the elsA software. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA, 13–17 January 2014; p. 0568. [Google Scholar]
  21. Othmer, C. A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows. Int. J. Numer. Methods Fluids 2008, 58, 861–877. [Google Scholar] [CrossRef]
  22. Othmer, C. Adjoint methods for car aerodynamics. J. Math. Ind. 2014, 4, 1–23. [Google Scholar] [CrossRef] [Green Version]
  23. Nemec, M.; Zingg, D.W.; Pulliam, T.H. Multipoint and multi-objective aerodynamic shape optimization. AIAA J. 2004, 42, 1057–1065. [Google Scholar] [CrossRef]
  24. Schramm, M.; Stoevesandt, B.; Peinke, J. Simulation and Optimization of an Airfoil with Leading Edge Slat; IOP Publishing: Bristol, UK, 2016; Volume 753, p. 022052. [Google Scholar]
  25. Rashad, R.; Zingg, D.W. Aerodynamic shape optimization for natural laminar flow using a discrete-adjoint approach. AIAA J. 2016, 54, 3321–3337. [Google Scholar] [CrossRef]
  26. Elham, A.; van Tooren, M.J.L. Discrete adjoint aerodynamic shape optimization using symbolic analysis with OpenFEMflow. Struct. Multidiscip. Optim. 2021, 63, 2531–2551. [Google Scholar] [CrossRef]
  27. Liou, M.S. A sequel to AUSM: AUSM+. J. Comput. Phys. 1996, 129, 364–382. [Google Scholar] [CrossRef]
  28. Liou, M.S.; Steffen Jr, C.J. A new flux splitting scheme. J. Comput. Phys. 1993, 107, 23–39. [Google Scholar] [CrossRef] [Green Version]
  29. Roache, P.J. Verification and Validation in Computational Science and Engineering; Hermosa: Albuquerque, NM, USA, 1998; Volume 895. [Google Scholar]
Figure 1. The algorithmic flowchart for the shape optimization procedure involves solving the flow field equations, formulating and solving the adjoint equations, computing the converged adjoint variables, and determining the optimization direction based on the adjoint solver variables.
Figure 1. The algorithmic flowchart for the shape optimization procedure involves solving the flow field equations, formulating and solving the adjoint equations, computing the converged adjoint variables, and determining the optimization direction based on the adjoint solver variables.
Applsci 13 04269 g001
Figure 2. (a) The Mach lines around the NACA 0012 airfoil under nonlinear flow conditions at M = 0.8 and α = 1 are presented. (b) Nonlinear pressure coefficient lines are shown around the same airfoil under the same flow conditions. (c) The pressure distribution on the surface of the NACA 0012 airfoil is presented for non-viscous flow at M = 0.8 and α = 1 , and the results are compared to the numerical results obtained by Liou and Steffen Jr. [28].
Figure 2. (a) The Mach lines around the NACA 0012 airfoil under nonlinear flow conditions at M = 0.8 and α = 1 are presented. (b) Nonlinear pressure coefficient lines are shown around the same airfoil under the same flow conditions. (c) The pressure distribution on the surface of the NACA 0012 airfoil is presented for non-viscous flow at M = 0.8 and α = 1 , and the results are compared to the numerical results obtained by Liou and Steffen Jr. [28].
Applsci 13 04269 g002
Figure 3. The residual of adjoint variables was calculated for 10 runs under the same mesh conditions, using the same input data as the example in the paper. These calculations were performed within a single loop consisting of 2400 iterations.
Figure 3. The residual of adjoint variables was calculated for 10 runs under the same mesh conditions, using the same input data as the example in the paper. These calculations were performed within a single loop consisting of 2400 iterations.
Applsci 13 04269 g003
Figure 4. The residual of adjoint variables was calculated for three different mesh conditions using the same input data as the example presented in the Multipurpose Optimization section. A single loop consisting of 2400 iterations was used for the calculation.
Figure 4. The residual of adjoint variables was calculated for three different mesh conditions using the same input data as the example presented in the Multipurpose Optimization section. A single loop consisting of 2400 iterations was used for the calculation.
Applsci 13 04269 g004
Figure 5. The figures illustrate various aspects of the optimization process of the multipurpose lift-to-drag ratio algorithm. (a) Changes in the cost function during the optimization process. (b) Represents the lift-to-drag ratio C l / C d obtained after applying the algorithm. (c) Shows variations in the pressure drag coefficient C d p , the viscous drag coefficient C d v , and the total drag coefficient C d during the optimization process. Finally, (d) displays the changes in the lift coefficient C l during the optimization process.
Figure 5. The figures illustrate various aspects of the optimization process of the multipurpose lift-to-drag ratio algorithm. (a) Changes in the cost function during the optimization process. (b) Represents the lift-to-drag ratio C l / C d obtained after applying the algorithm. (c) Shows variations in the pressure drag coefficient C d p , the viscous drag coefficient C d v , and the total drag coefficient C d during the optimization process. Finally, (d) displays the changes in the lift coefficient C l during the optimization process.
Applsci 13 04269 g005
Figure 6. The variations in the airfoil design parameters when applying the multipurpose lift-to-ratio optimization algorithm are shown, including the maximum thickness t, the location of the bend m, and the maximum bend p.
Figure 6. The variations in the airfoil design parameters when applying the multipurpose lift-to-ratio optimization algorithm are shown, including the maximum thickness t, the location of the bend m, and the maximum bend p.
Applsci 13 04269 g006
Figure 7. The proposed multipurpose optimization algorithm generated modifications in the pressure coefficient C p around the airfoil.
Figure 7. The proposed multipurpose optimization algorithm generated modifications in the pressure coefficient C p around the airfoil.
Applsci 13 04269 g007
Figure 8. The proposed multipurpose optimization algorithm generated modifications in the geometry of the airfoil curve.
Figure 8. The proposed multipurpose optimization algorithm generated modifications in the geometry of the airfoil curve.
Applsci 13 04269 g008
Figure 9. The Mach number contour around an airfoil was examined under the initial and final conditions of the proposed multipurpose optimization algorithm.
Figure 9. The Mach number contour around an airfoil was examined under the initial and final conditions of the proposed multipurpose optimization algorithm.
Applsci 13 04269 g009
Figure 10. The streamlines around the airfoil were observed in both the initial and final conditions of the multipurpose optimization.
Figure 10. The streamlines around the airfoil were observed in both the initial and final conditions of the multipurpose optimization.
Applsci 13 04269 g010
Table 1. The range of changes allowed in the determining components of the flap geometry during the optimization process.
Table 1. The range of changes allowed in the determining components of the flap geometry during the optimization process.
ParameterSymbolRatio
Airfoil Thicknesst0.1 c ∼ 0.4 c
Airfoil Maximum Camberm0%∼9.5%
Airfoil Maximum Camber Positionp0.1∼0.9
Table 2. The range of changes permitted in the determination of the components of the flap geometry in the optimization path.
Table 2. The range of changes permitted in the determination of the components of the flap geometry in the optimization path.
ParameterSymbolRatio
Airfoil Thicknesst0.1 c∼0.4 c
Airfoil Maximum Camberm0%∼9.5%
Airfoil Maximum Camber Positionp0.1∼0.9
Table 3. The flow conditions utilized in the proposed multipurpose shape optimization for the NACA4418 airfoil.
Table 3. The flow conditions utilized in the proposed multipurpose shape optimization for the NACA4418 airfoil.
ParameterSymbolRatio
Mach NumberM0.72
Reynolds NumberRe3000
Angle of Attack α 2.8
Target Lift Coefficient C l max 0.45
Target Drag Coefficient C d min 0.12
Table 4. The initial and final values of the components that determine the airfoil geometry and aerodynamic coefficients were optimized to achieve maximum lift and minimum drag coefficients.
Table 4. The initial and final values of the components that determine the airfoil geometry and aerodynamic coefficients were optimized to achieve maximum lift and minimum drag coefficients.
ParameterSymbolInitial ValueFinal Value
Maximum Thicknesst0.180 c0.143 c
Maximum Camberm4.000%3.986%
Maximum Camber positionp0.4000.397
Lift Coefficient C l 0.1010.398
Pressure Drag Coefficient C d p 0.1500.083
Viscous Drag Coefficient C d v 0.0300.037
Total Drag Coefficient C d 0.1790.120
Lift-to-Drag Ratio C l / C d 0.5643.328
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Tanabi, N.; Silva, A.M., Jr.; Pessoa, M.A.O.; Tsuzuki, M.S.G. Robust Algorithm Software for NACA 4-Digit Airfoil Shape Optimization Using the Adjoint Method. Appl. Sci. 2023, 13, 4269. https://doi.org/10.3390/app13074269

AMA Style

Tanabi N, Silva AM Jr., Pessoa MAO, Tsuzuki MSG. Robust Algorithm Software for NACA 4-Digit Airfoil Shape Optimization Using the Adjoint Method. Applied Sciences. 2023; 13(7):4269. https://doi.org/10.3390/app13074269

Chicago/Turabian Style

Tanabi, Naser, Agesinaldo Matos Silva, Jr., Marcosiris Amorim Oliveira Pessoa, and Marcos Sales Guerra Tsuzuki. 2023. "Robust Algorithm Software for NACA 4-Digit Airfoil Shape Optimization Using the Adjoint Method" Applied Sciences 13, no. 7: 4269. https://doi.org/10.3390/app13074269

APA Style

Tanabi, N., Silva, A. M., Jr., Pessoa, M. A. O., & Tsuzuki, M. S. G. (2023). Robust Algorithm Software for NACA 4-Digit Airfoil Shape Optimization Using the Adjoint Method. Applied Sciences, 13(7), 4269. https://doi.org/10.3390/app13074269

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop