Aerodynamic Shape Optimization of a NACA0018 Airfoil Using Adjoint Method and Gradient-Based Optimizer

. The purpose of this study is to demonstrate the suitability and efficacy of the adjoint method on the aerodynamic shape optimization on a simple symmetrical airfoil NACA0018 at low Reynolds number for wind turbine application. The adjoint method has been used in many pressure-based numerical simulations with various degrees of success leading to optimized geometries in their respective uses. ANSYS Fluent code was used in this simulation. Lift to drag ratio was defined as the observables for which adjoint sensitivities were formulated. The objective function of the optimization was set to maximize the lift to drag ratio of the airfoil by 20%. The optimization regime showed significant increase in lift and drag ratio from the initial baseline NACA0018 value of 0.0211 up to 3.66 for the optimal NACAOpt. The results demonstrate the potential of the adjoint solver paired together with the gradient-based optimizer to improve the geometry for shape optimization in many CFD applications.


Introduction
The vertical axis wind turbine (VAWT) is a mechanical device extracts wind energy. A typical Darrieus H-type VAWT has aerodynamically profiled blades that rotates along its vertical axis and is supported and connected to the main rotor shaft or central column by support arms. The main rotor shaft is then connected to the gearbox and generator located at the bottom of the VAWT and energy generated is transmitted to the grid for industrial or urban use. VAWTs utilizes airfoils to capture wind energy to spin the rotor and convert it to usable energy. Hence, the aerodynamics of a turbine blade plays an important role during turbine operation which generates research interest to maximize its efficiency.
The application of the adjoint method optimization procedure can be applied to the VAWT components to maximize efficiency. The adjoint method is an optimization procedure based on design sensitivities that has been used in various computational fluid dynamics (CFD) applications. Multiple studies have been presented showing its uses for shape optimization in various fields of studies and applications. The adjoint method encompasses calculations of the adjoint equations and residuals to obtain the design sensitivities for mesh morphing of the geometry which is then optimized using the tried and tested gradient-based optimization algorithm. In contrast to many other available optimization schemes, the adjoint method is able to handle a large number of design variables in a single adjoint calculation [1]. This means significant savings towards computation time and resources as well as better optimization efficiency with respect to the number of design variables.

Literature Review
With respect to aerodynamics and airfoil geometries, the H-Darrieus turbine blade can take the form of various different airfoil profiles. A turbine blade can take the shape of any symmetrical or unsymmetrical blade which each having its own advantages [2][3][4]. Many studies have been done on the VAWT to optimize or maximize its performance way by suggesting and comparing various different types of turbine blades. Comparisons between support arm geometries such as implementing symmetrical NACA blades [5,6] or variations of it [7,8] as well as other shapes [9,10] which generally tends to favour airfoils shapes when it comes to power efficiency losses. However, implementation of shape optimization techniques to generate a novel blade design for VAWT use has recently seen a rise in interest.
Optimization techniques with respect to the aerodynamic performance of a VAWT has been applied in various forms. Marinic-Kragic et. al. [11][12][13] has provided numerous examples of VAWT optimization performed on the Savonius type VAWT. An initial proposal for an optimization scheme using shape parameterization along with a genetic algorithm code in a CFD analysis of the VAWT [13]. The authors proposed a numerical framework to generate an optimized VAWT which was able to generate a new shape based on the input of the user and not based on an initial starting shape. A followup study by the same authors proposed a novel flexible shaped Savonius turbine that is able to obtain an increase of 8% in power efficiency [11] and later further examples of Savonius turbine optimization with a 12% improvement in a 2 bladed design compared to the baseline and close to 50% improvement for a 4 bladed and 6 bladed setup [12]. Elsaka et. al. [14] used response surface optimization together with the multi-objective Genetic Algorithm (MOGA) in a CFD simulation to find optimal turbine solidity of 0.31 and blade pitch angle of 3.6. The 2D results show a minimal improvement of about 2.9% over the baseline turbine but when comparing 3D simulations, the added parameter of aspect ratio resulted in a significantly higher percentage improvement of about 34.5% with the winglet design increasing turbine efficiency even higher especially at lower AR values. The applications of the adjoint method in geometry optimizations have been reported in many different applications. The adjoint based aerodynamic optimization of airfoils has been used for large scale transonic inviscid flow conditions primarily in aviation applications to reduce drag and has yielded good results [15][16][17]. Giles and Pierce [18] presented a simple adjoint method approach by demonstrating the optimization of an aeroplane showing significant reduction in drag coefficients. Mohebbi et. al. [19] showed the difference between a short iteration number and a longer one with the former run over 5 iterations obtained a 36.95% increase in lift to drag ratio while the later run over 20 iterations showed significantly higher increase by 93.34%. Fleichli et. al. [20] compared two optimization methods and found the Evolutionary Algorithm optimization achieved an average difference of 4.75% reduction between the inlet and outlet tube. The unrestricted adjoint solver however achieved a significant 19.23% improvement over the baseline. This is a difference of 14.48% between the two optimization methods. In the automotive industry, the adjoint method has been used for pressurebased optimization for various instruments. Horova et. al. [21] utilized the adjoint method to minimize the pressure loss of an intercooler filling line and achieved a significant 23.9% reduction and Lang [22] [10] illustrated the adjoint method optimization procedure and demonstrated the use of the adjoint solver to optimize various sections of a car related to its aerodynamics characteristics. In wind turbine optimization, multiple different studies have employed different optimization techniques. Dhert et. al. [23] implemented the on a HAWT turbine blade and obtained a torque increase of up to 22.4%. In VAWT, Day et. al. [24] employed the adjoint method optimization on a vertical axis wind turbine. The objective function was conservatively set to a 3% increase applied to the turbine blade tangential blade force coefficient and conducted at a relatively high TSR of 4. Improvement showed significant increase in tangential forces which resulted in higher lift and better efficiency. However, the general flow characteristic of the suggested turbine and performance curves such as power coefficient were not reported.
This study attempts to demonstrate the adjoint optimization method on a symmetrical NACA0018 airfoil at low velocity typically experienced by VAWT (Re ~1000). The objective function of the optimization is to increase the airfoil lift to drag ratio of the symmetrical NACA0018 and extract a novel shape. The result of this study is expected to be applied to the vertical wind turbine blade and support arm in future research.

Adjoint optimization
Aerodynamic shape optimization of the turbine support arm was done using the adjoint solver together with the gradient-based optimizer both included in the commercial CFD software ANSYS Fluent. An initial CFD calculation is conducted first with the operating conditions and inputs such as boundary conditions, initial geometry etc. which are the baseline for optimization. This process produces the initial baseline data such as XY plots of geometry lift and drag as well as flow field data such as velocity and pressure fields. Next, the observables, which are the lift, drag and lift to drag ratio, are then defined, and calculated as objective functions for improvements. For this procedure, the objective function was set to an increase of 20% to the lift to drag ratio. The adjoint solver was then run to gather sensitivity data of the mesh with respect to the observables and the initial flow field conditions. The sensitivity analysis is crucial as this process identifies regions of the mesh most sensitive to the defined observables and therefore the most suitable locations for shape deformation through node displacements in the mesh. Figure 2 shows the regions at the leading and trailing edge of the airfoil where mesh morphing is most sensitive to the changes in lift to drag. The mesh morpher was then initiated and, in the design tool tab, the bounding region defined to envelop the airfoil geometry. For the control points, 20 nodes were set for both the X and Y direction as shown in Figure 3. The gradient-based optimizer was initialized with a single condition of improving the lift to drag ratio by percentage increase of 20%. The mesh morphing and optimization was then started and conducted unrestricted up to 200 iterations until a suitable shape that meets the objective function was reached. A summary of the adjoint approach is presented in Figure  4. This study implements and compares two different objective functions. The first is to optimize the lift to drag ratio of the airfoil and the second to maximize the moment. The shape optimization is set to maximizing the objective function subject to the control points c. This can be summarized as The objective function J is defined as Equation (1) The value chosen for Langragian multiplier ϕ is such that it reduces the following equation to ∂ℐ ∂q = ϕ T ∂R ∂q Equation (2) And finally, the adjoint sensitivity equation is defined as,

Numerical simulation settings and computational domain
In this study, the governing equations were numerically computed using the ANSYS Fluent CFD software. Two dimensional, incompressible Reynolds-averaged Navier Stokes (RANS) approach was employed coupled with the Transitional SST turbulence model. Unstructured 2D rectangular mesh was used as shown in Figure 5 whereby the boundary layer was structured mesh. The velocity inlet was set to the curved edge on the left side of the domain and the pressure outlet set to the vertical edge on the right side. The remaining horizontal sides were set to standard wall with no slip conditions.  Mesh independency study was conducted to determine the domain size. As shown in Figure 6, it was decided that 9 x 10 5 mesh element size to be considered as the optimal balance between result accuracy and simulation time. This element mesh size was used for both NACA0018 and NACAOpt simulations. Second order discretization was selected for the pressure and momentum formulation and SIMPLE scheme was used for the velocity-pressure coupling. Table 1 summarizes the numerical parameters of the simulation.

Validation
Mesh validation was conducted to ensure the validity of the meshing procedure used in this study. The validation was conducted at Re 150,000 against results by Michna and Rogowski [25]. The validation was performed in similar condition to the referenced study. Shown in Figure 7, it can be observed that the differences between the current numerical predictions agree well with the referenced results. The trend is consistent between the current results and the referenced studies. Although the numerical optimization in this study adopts different airfoil conditions, the validation conducted above shows the meshing strategy employed has been validated to be as accurate as possible.  Figure 8 shows the geometric comparison between the NACA0018 and NACAOpt. It can be observed that the overall optimized NACAOpt geometry shows an unsymmetrical profile mainly showing most morphing on the lower side of the airfoil which was expected. The leading-edge profile scoped inwards and upwards at about 0.2 of the chord length. This results in a sharper leading-edge profile. The trailing edge lower side also shows an inward and upward morphing. Lastly, the overall profile of the NACAOpt was noticeably smaller than the NACA0018 airfoil.

Aerodynamic performance coefficients
The gradient-based optimizer was allowed to run until the it achieved convergence irrespective of the set goal increase of the objective function following a suggestion by [11]. This allows a significantly unrestricted calculation to obtain an improved performance. Table 2 compares the aerodynamic coefficients between the NACA0018 and the NACAOpt airfoil geometry at 0°. The lift coefficient and lift to drag ratio both experienced significant improvements of 191.3% and 190.44% increase respectively. The drag coefficient saw negligible changes. Figure 9 shows the comparison between the NACA0018 and the NACAOpt geometries. The NACAOpt shows significantly higher lift to drag ratio compared to the NACA0018 geometry at lower angles of attack between 0° to 12° degrees. Between 16° to 20° angle of attack, the differences between the two airfoils were negligible. This indicates that the NACAOpt generates better lift at lower angles of attack due to the optimized leading edge which promotes better tangential force acting upwards in the lift direction but then suffers from increased drag as the angle of attack increases. Figure 8. Comparisons of the initial NACA0018 geometry and the optimized NACAOpt geometry.  Figure 10 and Figure 11 shows the velocity and pressure contour comparison between the NACA0018 and NACAOpt. When comparing the velocity flow contour, while the stagnation point on the NACA0018 airfoil shown in Figure 10(a) was at the center of the leading edge, the stagnation point on the NACAOpt was found to have shifted towards the lower side of the airfoil shown in Figure 10(b). Flow separation happens towards the upper surface of the leading edge unlike the NACA0018 which happens at the centre of the leading edge. With the upper surface of the more streamlined and the lower surface of the airfoil cambered, the fluid travels faster over the upper surface while the cambered outline of the lower surface slows the fluid flow which delays flow reattachment. This phenomenon encourages a pressure difference between upper and lower surface which is favourable to lift.

Flow field contour comparisons
The NACA0018 pressure contour in Figure 11(a) predictably shows an even pressure distribution along its upper and lower surface whereas the NACAOpt in Figure 11(b) shows a much more chaotic and uneven distribution. The leading-edge experiences strong compressive pressure on the lower side with the upper surface of the blade experiencing suction similarly seen distributed along the lower surface of the airfoil. The symmetrical nature of the pressure distribution on the upper and lower surface of the NACA0018 would indicate a streamlined flow where the fluid flow would detach from the leading edge and reattach at the trailing edge at the same time hence contributing to poor lift at 0° angle of attack. Conversely, the uneven pressure distribution on the upper and lower surface of the airfoil due to the unsymmetrical shape of the NACAOpt. The leading edge experiences higher pressure at the leading lower edge surface would therefore contribute to its higher lift even at 0° angle of attack. The overall difference in magnitude of this pressure difference was found to be higher in the NACAOpt than the symmetrical NACA0018 airfoil.  Figure 12 shows the surface pressure coefficient comparisons along the upper and lower surface of the airfoils. From Figure 12(a), the NACA0018 airfoil sees a reduction in pressure coefficient up to about 0.1 of the x/c. before steadily increasing along the airfoil surface towards the trailing edge. The NACAOpt airfoil experiences a larger drop in pressure coefficient at the leading edge of the airfoil, which then increases at 0.1 x/c before steadily increasing towards the trailing edge like the NACA0018 airfoil. The overall pressure on the upper surface of the NACAOpt airfoil was noticeably lower than the NACA0018 airfoil. In Figure 12(b), the NACA0018 airfoil predictably shows a similar shape to that at the upper surface which is due to its symmetrical shape but the NACAOpt airfoil shows a more chaotic shape oscillating between positive and negative pressure coefficient along the lower surface of the airfoil towards the trailing edge. The pressure coefficients show a higher positive value between 0 to 0.2 x/c and a significant peak at 0.5 x/c but with noticeably lower troughs at about 0.3 x/c and 0.6 x/c. These locations are at significant curves along the airfoil profile on the lower side and can be attributed to the irregular cambered shape of the airfoil which contributes to the increase in lift generated.

Conclusions
The shape optimization process using the powerful adjoint method of a NACA0018 airfoil has been described. The adjoint method shape optimization scheme drastically improved the lift to drag ratio by 190.44% as well as lift coefficient by 191.30% respectively. Drag was found to have changed negligibly. The velocity and pressure contours showed a shift of the stagnation point higher as well as an uneven pressure distribution along the surfaces of the airfoil in the NACAOpt geometry. These results demonstrate the viability of the adjoint method in increasing aerodynamic performance of the airfoil. Continuing on from this study, the adjoint method will be implemented to optimize other output parameters such as airfoil drag and pitching moment as well turbine moment and power efficiency. These optimization schemes will be done on the turbine blade and support arm geometry for improved turbine performance.