Interactive Airfoil Optimization Using Parsec Parametrization and Adjoint Method

Featured Application: The potential application of the work presented in this article lies in the preliminary airfoil and wing design and optimization. It is in this area that the strengths of this approach (negligible computation cost, robustness, . . .) can be maximally exploited while the weaknesses that hinder it can be minimized. The code is ready to use for engineering practice. Abstract: In the development of interactive aerodynamic optimization tools, the need to reduce the computational complexity of flow calculations has arisen. Computational complexity can be reduced by estimating the flow variables using machine learning, but that approach has a number of hindrances. Avoiding these hindrances through lowering the computational complexity by stating the assumptions of inviscid incompressible potential flow is the focus of this article. The assumptions used restrict the applicability of this approach to only specific cases, but in engineering practice, these cases are quite widespread. The assumptions allowed the coupling of the adjoint method with parsec parametrization and the panel method, yielding a highly computationally efficient and robust tool for optimizing an airfoil’s lift coefficient ( C y ). The optimization of the NREL S809 airfoil was carried out, and the results were verified using the Xfoil 6.99 software. The Xfoil verification showed that by making minimal changes to the airfoil’s shape, the C y and lift-to-drag ratios were significantly improved. The improvement magnitude was over 94% for a 0 deg angle of attack (AoA) and over 16% for 6.2 deg AoA. This indicates an improvement in performance that is similar to that of some genetic algorithms, but with computational costs that are many orders of magnitude lower.


Introduction
The process used in the aerodynamic optimization of an airfoil relies on two main features.The first is a computationally efficient optimization method that can effectively minimize (or maximize) the cost function.The second feature, which is perhaps even more important than the first one, is the ability to represent an airfoil in terms of design variables (to parametrize the airfoil).Currently, there are many parametrizations for airfoils and many optimization methods, yielding dozens of possible combinations.However, the ability to perform CFD-based optimization reliably and reasonably quickly relies heavily on the discovery of the adjoint method in the 1970s and 1980s and its development in the 1990s and early 2000s.
Currently, the objective of aerodynamic optimization research and development is shifting towards fast interactive design optimization [1], where the results of the optimization have to be known almost instantly.The CFD-based adjoint optimization method, which is now an extremely well thought-out and proven concept that is implemented in both commercial (like Ansys Fluent) and open source (such as OpenFOAM) big CFD codes, is not well suited to this task.This is due to the vast computational complexity of modern-day CFD simulations, which makes the optimization too slow to be considered interactive, even with the use of the most efficient optimization methods [1].It follows that the computational complexity of both the flow solution and the optimization algorithm must be substantially reduced to allow interactive design optimization [1].This can be conducted in at least two ways.The first is presented in [1] and its references and involves the use of modern machine learning (ML) methods to estimate the flow variables (and sometimes even adjoint variables when coupled with the adjoint method) and force coefficients.The other way, which is presented in this article, is to stick to the adjoint method for optimization while making some assumptions about the flow, which allows the use of much simpler (and thus faster) methods for solving the flow variables, such as the panel method.This approach (when cautiously formulated) almost completely avoids three main hindrances of the ML approach: extrapolatory predictions, the curse of dimensionality, and the need for vast sets of good-quality training data [1].However, in the adjoint-based assumption approach, great care has to be taken when placing the design point and simplifying the flow, as some important flow phenomena may be omitted in the simplification.This means that for some problems, especially where compressible and/or viscous flow must be assumed, it may be more feasible to use ML methods to deal with the challenges of interactive optimization.
The adjoint method is a computationally efficient optimization method, and thanks to its efficiency and accuracy, it is well suited for optimization in fluid mechanics, as stated in [2,3].However, its complexity in both the mathematical apparatus and the programming implementation made initial progress after its proposal rather slow-paced [4].Since its proposal, it has been successfully used in a variety of optimization problems in the field of fluid mechanics (both external and internal), such as those in [2][3][4][5][6][7].The adjoint method comes from the theory of optimal control, as stated in [3,8], and was first proposed by Olivier Pironneau (*1945) and Anthony Jameson (*1934) [7].Jameson developed and successfully used the first optimization code based on the adjoint approach [9].The adjoint method, unlike the finite difference method, evaluates the gradient of the cost function indirectly using the adjoint variables [2, 3,8].Consequently, the computational effort tied to the computation of the gradient of the cost function is not coupled to the number of design variables and is always roughly equal to the computation of two flow solutions [2, 3,8].Despite the computational efficiency of the adjoint method, it is good practice to parametrize the optimized geometry with an appropriate number of design variables and to avoid the use of redundant ones, thereby preventing (or at least minimizing) problems related to the geometrical or physical infeasibility of the solution obtained.This is particularly applicable when the original problem has to be simplified for optimization.
To obtain the design variables required by the optimization method, there are, at least in the case of an airfoil, literally dozens of parametrization methods [10,11].From those, parsec parametrization is only one of the many well-proven parametrizations available.Those parametrizations involve 4-and 5-digit NACA series and their modifications [12], class shape transformation (CST) [13], the analytic equation, as presented in [14], and so on.From those parametrizations, the 4-and 5-digit NACA series and analytic equation from [14] have less than or equal to six parameters.CST can have an arbitrary number of parameters, and parsec has 11 or 12, depending on its variant [13,15].There are even some modified parsec methods [15].The parsec method has been extensively used for airfoil optimization, such as in [13,[16][17][18].It has strong control over the airfoil's leading edge, crest locations, and curvatures at the crest locations [15].On top of that, its parameters have a clear and straightforward geometrical meaning, giving better insight into the airfoil shape.Parsec is a very powerful parametrization for covering general airfoil shapes with a reasonable amount of design variables [10,19].Its performance in geometrical accuracy is roughly comparable to that of the Bezier, RBF, and Hicks-Henne parametrizations, with a similar number of design variables [11].In contrast, its performance in the accuracy of aerodynamical properties is comparable to all the other parametrization methods covered in [11], with a similar number of design variables.However, parsec lacks precise control of the trailing edge shape [10], but this can be beneficial in this case, as the trailing edge shape has a very influential position in the panel method that may not be related to real physical phenomena.
For this article, the panel method proposed by Hess and Smith in 1966 [20] was used as a flow solver, allowing a very low computational cost.The simplicity of the panel method means that great care must be taken when selecting the parametrization.The ability to prescribe a nearly arbitrary trailing edge geometry (which is possible with many parametrizations) may hinder the practical usability of the optimization, either through manufacturing difficulties or through the flow phenomena the panel method does not account for.On top of that, too many design variables or excessively strong control over the trailing edge shape can lead to geometrically impossible (distorted) shapes arising during optimization.The parsec parametrization was used because a proven parametrization method was needed, and regarding the flow solver used, the assurance of the geometrical and physical feasibility of the resulting shape was required.
The next thing to consider is the uncertainty of the input data, together with the finite manufacturing precision.The dependence of the optimal shape on the input data uncertainties, together with the influence that the finite manufacturing precision has on the performance of both the original and optimized airfoils would be a very interesting research topic.Many methods for assessing the influence of uncertain operating conditions on the solution exist.Those methods and their possible usage can be found for example in [21,22].This topic, however, is so broad that a full-length article would be needed to cover it in appropriate detail; thus, it will not be discussed here.
The main goal of this article is to study the possibility of creating an interactive airfoil optimization procedure that runs very quickly; this procedure is based on an adjoint method coupled with parsec parametrization.The second goal is to provide Xfoil data for verification of the optimization's efficiency.This article also focuses on comparing the efficiency of the created optimization code and its results with other articles that discuss 2D airfoil optimization.The optimization algorithm created in this research can be of great value in all fields that use airfoils.These include not only aircraft and wind turbine design and development but also vehicle dynamics and stability, as can be seen, for example, in [23].

General Basics of Optimization
The optimization methods vary in terms of their computational complexity, cost function requirements, type of search (local or global), etc.In fluid mechanics, great emphasis is placed on computational efficiency, as every solution of the flow field is very computationally expensive [1,3,9], and the computational costs of the optimization methods vary by many orders of magnitude.Global methods are, for this reason, rather sparsely used in fluid mechanics, and the far better choices are usually local methods, as they are far less computationally demanding [1].Despite their locality, lower computational costs make them more suitable for most CFD optimization problems, and the locality of the search is not a problem when the cost function is carefully defined [3].The computational complexity of various methods can be seen in Table 1, which underlines the computational efficiency of the adjoint method.The simplex method seems even more efficient, but Table 1 does not provide complete information about the computational demands of the simplex method.The data in Table 1 were obtained from [2] and as the results of a theoretical analysis of optimization algorithms.Those data are valid when the algorithms proceed normally.In the case of the simplex method, however, the startup of the method is computationally more expensive, and the computation cost of the first iteration is equal to m + 1 (simplex in m-dimensional space has m + 1 vertices).On top of that, something called "stalled convergence" can be expected to happen in the simplex method.Stalled convergence means that the algorithm gets stuck, possibly far from the desired extremum.Solving stalled convergence in the case of the simplex method means rescaling the simplex; thus, the computational complexity of evading the stalled convergence is equal to m (one vertex from the original simplex can be kept).This means that the computational complexity from Table 1 is the absolute best-case scenario for the simplex method, and the real computational complexity may be much higher.The other algorithms do not particularly suffer from this problem.The symbols used in Table 1 are defined as follows: m. . .number of parameters (dimensionality of design space); k. . .number of individuals in each population (usually between 10 and 1000).

Adjoint Method: Continuous and Discrete Approach
There are two main approaches to adjoint optimization: continuous and discrete [3,4].They differ in the sequence of mathematical operations used during their derivation.Each of those two approaches can then be treated from either a Lagrange or a duality viewpoint as both yield the same results [7].In this paper, the Lagrange viewpoint is used as it connects more naturally to the constrained optimization.
The continuous approach favored by Jameson seems more natural from the perspective of physics and continuum mechanics.In this approach, the cost function of a continuous problem is defined first.Then, the adjoint partial differential equations (PDEs) are derived from the continuous flow equations, together with their boundary conditions (BCs) and possible initial conditions (ICs).All the equations are then discretized separately and numerically solved.After solving those equations, it is possible to compute the gradient of the cost function with respect to the design variables.In this approach, great care has to be taken when deriving the adjoint equations, especially their BCs [3].This makes this approach very complex and challenging in terms of math, although it can be computationally more efficient than the discrete version and can have lower memory requirements [4].The continuous version is also critical for understanding the behavior of the adjoint equations and their physical significance and for assessing the possibility of ill-posed problems [3,4].The continuous adjoint method workflow is shown in Figure 1.The discrete adjoint method differs in terms of the order of operations.The beginning is the same as in the continuous approach, i.e., the definition of the cost function in the continuous problem.After that, the continuous problem is discretized as a whole.This means the discretization of the fluid flow PDEs with the corresponding BCs and ICs implemented and the formulation of the discrete analog to the continuous cost function as a function of many variables.Only then is the adjoint set of algebraic equations derived from the discretized flow PDEs.The advantages and disadvantages of this approach can The discrete adjoint method differs in terms of the order of operations.The beginning is the same as in the continuous approach, i.e., the definition of the cost function in the continuous problem.After that, the continuous problem is discretized as a whole.This means the discretization of the fluid flow PDEs with the corresponding BCs and ICs implemented and the formulation of the discrete analog to the continuous cost function as a function of many variables.Only then is the adjoint set of algebraic equations derived from the discretized flow PDEs.The advantages and disadvantages of this approach can be found in [4].This approach avoids much of the complexity of the continuous approach, but great care must be taken when discretizing the problem.The discrete adjoint method workflow is shown in Figure 2.

Discrete Adjoint Method from the Lagrange Point of View: Governing Equations
In this article, the nomenclature defined by Jameson in [3] is used.Thus,  is a set of dependent variables, which are usually velocity components at the integration points, pressure at the integration points, and the like.In the panel method, the meaning of  is a bit more complex and is clarified later. is the set of independent parameters (design variables) undergoing optimization; in this case, these are the parsec parameters describing the shape of an airfoil.Vectors  and  are implicitly tied via equations R, which allows the computing of  from the known parameters .The fixed parameters (like viscosity, etc.) are not mentioned here as they do not change during optimization.With all the variables defined, it is possible to define the discrete version of the cost function I, where  and  are subject to constraints R.This leads to the classical constrained extremum problem solved via the Lagrange multipliers .
(, , ) = (, ) +  • (, ) (2) From the gradient, the total differential of  can be defined.However, as all variables in the function  are treated as independent in the Lagrange approach, the gradient of  with respect to  can only be solved when the first and last terms in Equation ( 4) are equal to zero.From there, two systems of Equations ( 5) and ( 6) arise, containing 2n equations in total, which makes it possible to compute  and .

Discrete Adjoint Method from the Lagrange Point of View: Governing Equations
In this article, the nomenclature defined by Jameson in [3] is used.Thus, w is a set of dependent variables, which are usually velocity components at the integration points, pressure at the integration points, and the like.In the panel method, the meaning of w is a bit more complex and is clarified later.F is the set of independent parameters (design variables) undergoing optimization; in this case, these are the parsec parameters describing the shape of an airfoil.Vectors w and F are implicitly tied via equations R, which allows the computing of w from the known parameters F. The fixed parameters (like viscosity, etc.) are not mentioned here as they do not change during optimization.With all the variables defined, it is possible to define the discrete version of the cost function I, where w and F are subject to constraints R.This leads to the classical constrained extremum problem solved via the Lagrange multipliers λ.
From the gradient, the total differential of ϕ can be defined.However, as all variables in the function ϕ are treated as independent in the Lagrange approach, the gradient of ϕ with respect to F can only be solved when the first and last terms in Equation ( 4) are equal to zero.From there, two systems of Equations ( 5) and ( 6) arise, containing 2n equations in total, which makes it possible to compute w and λ.
In the case of linear constraints R, Equation ( 5) can be written in matrix form, and the Jacobi matrix on the left side of Equation ( 6) can be easily computed.The system of equations then forms two systems of linear algebraic equations (LAEs): ( 8) and (9).
After solving for w and λ, it is easy to obtain the gradient of ϕ with respect to F. It can be proven that this gradient is identically equal to the gradient of I, as the constraints R must always be satisfied (5).Then, the simple algorithm of the steepest descent is applied to obtain the new approximation of the optimal parameters (11).

Parsec Parametrization
Parsec parametrization was chosen for the optimization as it provides good accuracy for general airfoil shapes in terms of geometry and fluid flow characteristics [10,11].Its parameters also have clear and straightforward geometrical meaning, allowing easy fixation of certain geometric quantities and better engineering insight into the airfoil shape.In the basic parsec parametrization, there are 11 parameters [10].For the optimization, a slightly modified version of parsec, presented in [13], was used, allowing better control of the airfoil's shape near the leading edge [13].The parameters are very similar to the ones used by the basic parsec; the only difference is that the leading edge radius is prescribed for the upper and lower surfaces separately.The geometrical meaning of those parameters is shown in Table 2 and Figure 3.The shape of the airfoil is then described by Equations ( 12) and ( 13) as a linear combination of base functions [10].The coefficients for this linear combination are computed using Equations ( 14) and ( 15), as stated in [13].
Table 2. Parsec parameters and their meaning.

Basic Parsec Parametrization [10]
Parsec Variant Used for Optimization [13] Parameter Meaning Parameter Meaning ------r lo [1] Lower surface leading edge radius x lo [1] Lower surface crest location-horizontal position x lo [1] Lower surface crest location-horizontal position y lo [1] Lower surface crest location-vertical position y lo [1] Lower surface crest location-vertical position y xxlo [1] Lower surface crest curvature y xxlo [1] Lower surface crest curvature r le [1] Leading edge radius r up [1] Upper surface leading edge radius x up [1] Upper surface crest location-horizontal position x up [1] Upper surface crest location-horizontal position y up [1] Upper surface crest location-vertical position y up [1] Upper surface crest location-vertical position y xxup [1] Upper surface crest curvature y xxup [1] Upper surface crest curvature Trailing edge direction Trailing edge wedge angle Trailing edge wedge angle y te [1] Trailing edge location-vertical position y te [1] Trailing edge location-vertical position ∆y te [1] Trailing edge thickness ∆y te [1] Trailing edge thickness Trailing edge wedge angle  [1] Trailing edge location-vertical position  [1] Trailing edge location-vertical posi Δ [1] Trailing edge thickness Δ [1] Trailing edge thickness It has been assumed that the trailing edge thickness (∆y te ) is always equal to zero (sharp trailing edge) and is therefore not a design variable.This may limit the applicability in some cases, but the benefits in terms of the robustness of the algorithm were deemed to outweigh this drawback.The sharp trailing edge makes it easier to properly enforce the Kutta condition and prevents the optimization from making illegal changes to the trailing edge geometry.The area around the trailing edge is where the viscous effects are usually the most profound since the wake tends to start forming there.Because the optimization cannot account for those viscous effects (and especially the wake formation), it was decided to restrain it in this area as much as possible.On top of that, the trailing edge can be tweaked in many ways, most of them far exceeding the capabilities of the parsec parametrization, as the lack of trailing edge control is deemed to be one of its main weaknesses [10].For the trailing edge fine-tuning, methods capable of properly capturing viscous phenomena must be used.Vector F of the 11 design variables based on parsec parametrization is given as follows.F = r lo , x lo , y lo , y xxlo , r up , x up , y up , y xxup , α te , β te , y te T (16)

Potential Flow and Panel Method
For modeling the flow around an airfoil, the inviscid incompressible potential flow assumption was used.This assumption may be rather restrictive but allows the use of potential flow theory and the panel method described in [24].The incompressibility condition may be weakened in this case by the fact that the optimization presented in Section 4 holds, as long as the Prandtl-Glauert compressibility correction rule is valid.According to [25], this can be stated for flows with Ma ∞ ≤ 0.5, although great caution should be taken when dealing with flows with Ma ∞ > 0.4 in this manner.To obtain the potential flow solution, the Hess-Smith panel method [20] was used.An explanation of the governing equations of potential flow theory, together with the theory and governing equations of the panel method, can be found in [20,24].After employing the equations from [20,24], the lift coefficient (C y ) can be defined and rewritten in dimensionless variables (17); the definition of the dimensionless variables can be found in (18).
For the easier derivation of the equations in the following sections, the chord length L can be factored out from the integral in (17) and transferred to the left side of the equation, yielding the relationship for the continuous cost function (19).
Now, it is possible to discretize the whole problem, including the governing harmonic equation for velocity potential, as presented in [24] and the cost function (19).The discrete version of the cost function is obtained simply by discretizing the integral in (19).To obtain the discretized flow equations, a Hess-Smith panel method [20] is used.The shape of the whole airfoil is discretized into many line segments (panels), leading to the panel method equations, as derived by [20].The dependent variables w in the optimization are the source on each panel and the circulation around the airfoil.The governing equations of the panel method are presented below, together with the discrete version of the cost function.The subscript i is the index of the panel.The panel numbering is shown in Figure 4.
Aw = b, w = [q 1 , . . . ,q N , Γ] T (22) Appl.Sci.2024, 14, x FOR PEER REVIEW The exact definition of the elements of matrix A and vector b can be found From that definition, it can be seen that the value  can be factored out of each the vector b, and the whole system of LAEs (22) can be divided by it, yielding Eq (23).
The tangent velocity on each panel can be extracted directly using the panel m [20].That relationship can then be divided by  , yielding Equation (24), which substituted into the cost function (20).This last modification completes the arran of equations needed to derive the adjoint optimization.

Adjoint-Based Optimization with Panel Method and Parsec Parametrization
To obtain the equations for the adjoint optimization, it is necessary to start with tions ( 9) and (23).Matrix A and vector b are known from [20].The only part that m specified in Equation ( 9) is its right-hand side, which contains, in this context, the g of I with respect to  .The exact definition of the elements of matrix A and vector b can be found in [20].From that definition, it can be seen that the value c ∞ can be factored out of each term in the vector b, and the whole system of LAEs ( 22) can be divided by it, yielding Equation (23).
The tangent velocity on each panel can be extracted directly using the panel method [20].That relationship can then be divided by c ∞ , yielding Equation (24), which can be substituted into the cost function (20).This last modification completes the arrangement of equations needed to derive the adjoint optimization.

Adjoint-Based Optimization with Panel Method and Parsec Parametrization
To obtain the equations for the adjoint optimization, it is necessary to start with Equations ( 9) and ( 23).Matrix A and vector b are known from [20].The only part that must be specified in Equation ( 9) is its right-hand side, which contains, in this context, the gradient of I with respect to ŵ.
In (25), C is defined as follows: The definition of the gradient comes directly from (10), and its constituent derivatives are specified in (28) and (29).
Now, it is possible to continue processing the individual terms in Equations ( 28) and ( 29).The derivatives of the vectors ĉt and b, and matrix A can be obtained semianalytically, greatly reducing the number of numerical derivatives computed, which saves a lot of computational costs.It can be shown that only the following derivatives must be computed numerically: To evaluate these derivatives, the first forward difference is used.This approximation is first-order accurate but computationally very cheap.Because the code runs in MATLAB R2023b, which uses double precision as standard, a rather small step when evaluating the finite difference can be used without losing precision to rounding errors.It has been found that the value δ = 10 −8 gives accurate results and a stable optimization algorithm.The value δ = 10 −9 has been tried out without any noticeable improvement in precision or stability.All the numerical derivatives are computed as presented in (31).
Now, it is possible to define the following auxiliary variables.These variables arise when obtaining the analytical derivatives of individual terms in ( 28) and (29).They are not essential for the math, but they make programming implementation easier and far cleaner.
With the use of the variables defined in (32)-(36), it is now possible to rewrite all the terms in Equations ( 28) and (29) in the following manner: Then, the gradient of the cost function with respect to F can be constructed, and the new approximation of the optimal design variables F can be computed.In the gradient G, ŵ is the solution of ( 23), and λ is the solution of (9), where the gradient term on the right-hand side is given in (25).
The value ζ is the iteration step.For stability reasons, it proved beneficial to define it in the manner shown in (46).This definition is equal to using a unit gradient, which eliminates the stability problems tied to the strong nonlinearity of the cost function in some regions of the design space.The value ζ 0 is the nominal iteration step and is user-defined.

Resulting Shape Feasibility and Convergence of the Optimization
The programming implementation of the procedure presented in Section 4.1.and Section 4.2. is rather simple but requires the addressing of a few tasks.The first task is to prevent the occurrence of distorted shapes during optimization.In the solution of this problem, the clear geometrical meaning of parsec parameters is very handy.With that meaning, it is possible to simply prescribe the conditions that assure the geometrical feasibility of those shapes.This is carried out by monitoring the values of the parameters (like r up , r lo , x up , x lo , and β te ); if they get outside the allowed range, the optimization is immediately stopped.The allowed range can be completely controlled by the user, enabling custom restrictions to be made.On top of that, the airfoil paneling is constructed in such a way that the panel nodes in the top and bottom halves of the airfoil share x coordinates.Because of that, it is easy to compare the y coordinates of the corresponding nodes and to make sure that the y coordinate of the upper-surface node is always greater than that of the lower-surface one.If any self-intersection of the airfoil's surface is detected in this way, the optimization is immediately stopped, providing another layer of protection besides the allowed parameter values.
The next task is convergence monitoring.In the case of this optimization, monotone convergence is required.In the case of the worsening of the value of the cost function between subsequent iterations, a warning is displayed, but the optimization process is not terminated.This allows the procedure the possibility to regain monotone convergence once again.
The final task is to define the convergence criteria.As the inviscid incompressible potential flow is assumed, the C y may become unbounded, as there is no flow separation allowed other than that at the trailing edge.This means that reaching the true maximum of C y in the incompressible inviscid scenario may be, and probably is, impossible.This renders the classical convergence criteria for reaching the extremum unusable and also means that geometrically feasible but nonsensical and obviously wrong shapes are possible if the optimization is allowed to continue for too many iterations (more than a few tens).To stop this from happening, other convergence criteria were implemented.The optimization has three stopping criteria in total.The first states that convergence is reached when a defined improvement, whether absolute or relative, has been made.The second one stops the iteration when a defined change in shape, measured as the root mean square difference from the original shape, has been made.The final one limits the maximum number of iterations allowed.
The proposed criteria, however, cause the optimal result to be influenced by the threshold values chosen.The shape change criterion is the most objective of the three, as it can be interpreted as defining the region of design space for the algorithm to search in.It can also be tackled from the engineering viewpoint, as the baseline shape is usually chosen for a reason, and one may not want to change it much during the optimization.For more objective criteria to be defined, it would be possible to add a simple boundary layer solver based on the von Karman momentum integral equation and Pohlhausen's method to determine whether the flow separation was occurring on the airfoil.The stopping criterion could then be the occurrence of the separation, for example (or, as a variation of this criterion, the reaching of a certain distance from the trailing edge).This modification, however, would bring with it some computation costs, which in this context would be rather large.
The script conducting the optimization proposed in this paper was written in MATLAB R2023b.The script operates according to the flowchart shown in Figure 5.

Adjoint Optimization Results: Incompressible and Inviscid Flow
The most important results of the optimization, using the example of the NREL S809 airfoil, are presented below in Table 3 and Figure 6.Different airfoils were tried, yielding very similar results.These results show that massive improvement in  can be obtained by minimal change in the design variables, which implies minimal changes in the geometry.These minimal changes to the shape allow the assumption that the drag coefficient

Adjoint Optimization Results: Incompressible and Inviscid Flow
The most important results of the optimization, using the example of the NREL S809 airfoil, are presented below in Table 3 and Figure 6.Different airfoils were tried, yielding very similar results.These results show that massive improvement in C y can be obtained by minimal change in the design variables, which implies minimal changes in the geometry.These minimal changes to the shape allow the assumption that the drag coefficient will not be affected.This was later proven true by the Xfoil 6.99 computation results.The changes in the design variables and C y from the optimization solver for the nominal angle of attack (AoA) of 0 and 10 deg are presented in Table 3. Optimization settings with 50 iterations and a nominal iteration step ζ 0 = 0.0002 (a very conservative step choice) were used throughout this study.The results have been rounded to four significant digits in the case of the design variables and four decimal places in the case of C y .
Table 3. Parsec parameters of NREL S809 airfoil before and after the optimization.

Xfoil Verification Computations: Compressible and Viscous Flow
The results presented in this section were obtained using the well-established Xfoil 6.99 software.The Reynolds and Mach numbers were set to Re = 750,000 and Ma = 0.02, allowing a direct comparison between the obtained results and those published in [13].The nominal AoA for the adjoint optimization was zero; the resulting values of the design variables can be found in Table 3, Col. 4, and the original and optimized shapes can be seen in Figure 6.The dependency of the optimal shape on AoA is, at least in this case, very weak, as can be seen in Table 3.This means that the value of AoA is not of great importance in this case and that it has negligible effects on the optimization result.The result's lack of dependence on AoA, however, cannot be generalized and has to be assessed on each optimized airfoil separately.
The adjoint optimization results presented in Table 4 were compared to those acquired by Akram and Kim [13], as shown in Table 5.For details regarding the acquisition of the data in Table 5 please see the original paper by Akram and Kim [13].The CFD computation model, grid, BCs, and other necessary information used for generating the CFD data in Table 5 are described in great detail in [13].A comparison of the different computation methods used during the research presented in this paper can be seen in Table 6 which directly compares those simulation methods to the available experimental data.For the viscous Xfoil 6.99 simulation, Re = 750,000 and Ma = 0.02 were used to match the experimental data.From this comparison, it can be concluded that the inviscid computation methods consistently overestimate the lift coefficient, which is an expected aspect of their behavior.Also, the accuracy of the inviscid methods can be considered sufficient only when there is no large flow separation occurring since they are unable to predict and account for this phenomenon.If a large separation occurs, their results can qualitatively differ from reality.Those methods cannot be used for computing the drag coefficient; thus, they cannot compute the L/D ratio either.However, as their results are consistent, the MATLAB (R2023b) version can be used as the highly efficient core for adjoint optimization, as long as no large flow separation is expected to occur.It can generally be said that a significant separation starts occurring when the lift coefficient curve starts deviating from the linear progression seen near zero AoA (roughly 6-8 deg AoA for NREL S809 airfoil at the presented values of Re and Ma).As for the compressible and viscous simulation, it can be concluded that this Xfoil 6.99 computation is accurate enough to serve as a verification computation for the optimization in this case.Thus, the benefit of the optimization can be assessed by the compressible and viscous simulation in Xfoil 6.99.When assessed in this way, the benefits should be maintained in the real world.such as those in [26,27].The experimental data presented in Figures 7 and 8 were taken from [26].The computations were carried out using the compressible and viscous Xfoil 6.99 simulation.

Conclusions
This article focused on the fast interactive airfoil shape optimization achieved by assuming an incompressible inviscid potential flow.The objective of the optimization was to improve the lift coefficient of the airfoil.The optimization was carried out on the well-known NREL S809 airfoil, which is used on wind turbines.From the results obtained, it can be concluded that coupling the adjoint method with parsec parametrization and the panel method can yield feasible results and can run at a speed that is more than fast enough for it to be considered interactive.The panel method and the adjoint optimization code were written in MATLAB R2023b, and the verification computations were carried out in Xfoil 6.99.From the results obtained, the following statements can be derived: 1.
The optimization showed a significant improvement to the objective function C y .For 0 deg AoA, the improvement was 94.7%, and for 6.2 deg AoA, it was 16.1%.2.
A big improvement was also observed in the L/D ratio, which matched or even exceeded the improvement to C y , as the optimization did not noticeably affect the C d .3.
For different AoAs (0 and 10 deg), the resulting shape was almost the same, meaning that the shape optimized for one operating point should work well in a wide range of AoAs, adding to the practical usability of the optimized airfoil.4.
The runtime of the adjoint optimization was, for the same airfoil with the same cost function, orders of magnitude shorter than the runtime of the GA coupled with Xfoil while yielding similar results.
Overall, this coupling seems very promising, but it is not without its difficulties.For example, further work would be needed to incorporate additional constraints into the optimization, most likely via the penalty formulation.However, the approach presented in this paper appears to suit those engineering problems where the inviscid incompressible flow can be assumed to be better than the ML approach because it avoids most of the main hindrances of ML while yielding convincing results.

Figure 1 .
Figure 1.Flowchart of the continuous adjoint method.

Figure 2 .
Figure 2. Flowchart of the discrete adjoint method.

Figure 3 .Figure 3 .
Figure 3. Parsec parameters of an airfoil and their geometrical meaning.

Figure 4 .
Figure 4. Node and panel numbering in the panel method.

Figure 4 .
Figure 4. Node and panel numbering in the panel method.

Figure 5 .
Figure 5. Flowchart of the optimization script.

Table 1 .
Computation cost of 1 optimization iteration-compared to single flow solution [2].