This notebook demonstrates how to use a common method found in the literature of airfoil optimization, parametric optimization. To this end we try a gradient-based method from the scipy optimization toolbox and use virtual displacements to allow the rapid the evaluation of the different geometries.
Do not take this notebook as a best-practice notebook for gradient-based optimization - unless you for some reason need or prefer the constraint of a PARSEC parameter space. With gradient based optimization, there are better ways to avoid the pitfalls of airfoil optimization [1] than to severely restrict the design space using parametrization, e.g. constraining the resulting geometry to avoid high curvature. And there are more parametrizations available that offer a wider range of possible geometries.
Here, we restrict the design space to a PARSEC [2] parametrization and only look for symmetric airfoils, which reduces the parameter space to 5 values. We look at optimal airfoil shapes for drag at different Reynolds numbers and airfoil surface areas. This is a problem that can be solved by gradient-free optimizers as well as by gradient based solvers.
There are a lot of cases in the literature that are concerned with the PARSEC parameterization and optimizing airfoils. For example, in [4] the PARSEC parametrization is used together with XFOIL [5] to find optima under different constraints and weights using heuristic search methods.
Disclaimer: You see here a case where most of the points where we look for optima converge to a sensible optimum, and where it does not the final iterate is perfectly sensible, too. This will not be the case when you just enter some random things. You will need to think about your goal function, how to scale it, which algorithm to use, what you need in terms of viiflow tolerance and so on - and sometimes luck.
import os
os.environ["MKL_NUM_THREADS"] = "1"
import viiflow as vf
import viiflowtools.vf_tools as vft
import viiflowtools.vf_plots as vfp
import numpy as np
from numpy import pi
import matplotlib
import matplotlib.pyplot as plt
import logging
import time
logging.getLogger().setLevel(logging.ERROR)
%config InlineBackend.figure_format = 'svg'
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = [12, 6]
# The code to generate parsec airfoils has been copied from
# https://github.com/dqsis/parsec-airfoils
# by Dimitrios Kiousis. The code is part of the git repository of this notebook.
from parsec import parsec_airfoil
# Sample coefficients (NACA0012)
# The parameters are from [3]
rle = .015
x_pre = 0.29663
y_pre = -0.06002
d2ydx2_pre = 0.45150
th_pre = 10*pi/180
x_suc = 0.29663
y_suc = 0.06002
d2ydx2_suc = -0.45150
th_suc = -10*pi/180
yte = .0
N = 140
param = np.asarray([rle,x_pre,y_pre,d2ydx2_pre,th_pre,x_suc,y_suc,d2ydx2_suc,th_suc])
parsecfoil,_ = parsec_airfoil(N,yte, rle,
x_pre, y_pre, d2ydx2_pre, th_pre,
x_suc, y_suc, d2ydx2_suc, th_suc)
Base = vft.repanel(parsecfoil,N)
BaseAirfoil = Base
We not only need the PARSEC parameter definition of an airfoil but also the gradients of all quantities of interest with respect to the parameter.
Viiflow returns the paraemters with respect to the virtual displacements, that is the shift of the geometry normal to its current surface at every discretization point.
To calculate the gradient with respect to the PARSEC parameters, we use the viiflowtools function virtual_displacement_from_geometry(Base,Target), which calculates the displacement given a desired shape (Base) from a given shape (Target).
The gradient is then given by finite differences of this function.
For symmetric airfoils, the parameters needed are reduced to 5.
# Build airfoil from parameters and, if needed, calculate finite difference gradient of displacement-from-parameters
def dDeltadParam(param,BaseAirfoil,calc_gradient):
N = BaseAirfoil.shape[1]
Airfoil0, dAFy = parsec_airfoil(N,0,*param,gradients=calc_gradient)
[delta,Shape,delta_X,delta_Y] = vf.virtual_displacement_from_geometry(BaseAirfoil,Airfoil0)
dDeltadparam = None
if calc_gradient:
dDeltadparam = delta_Y@dAFy[:,1::] # Only parameters from yte onwards
return [delta,dDeltadparam,Shape,Airfoil0]
def GetNormals(X):
# Calculate Normals from neighbours
Normals = np.zeros(X.shape)
NB = Normals.shape[1]
Dif = np.diff(X)
DifC = Dif[0,:] + 1j*Dif[1,:]
NC = np.zeros((NB,))+1j
NC[0] = DifC[0]*1j
NC[1:NB-1] = (DifC[0:NB-2]/np.abs(DifC[0:NB-2])+DifC[1:NB-1]/np.abs(DifC[1:NB-1]))*1j/2
NC[NB-1] = DifC[NB-2]*1j
Normals[0,:] = -np.real(NC)/np.abs(NC)
Normals[1,:] = -np.imag(NC)/np.abs(NC)
S = np.cumsum(np.r_[0,np.abs(DifC)])
return [Normals,S]
# For symmetric airfoils
def dDeltadParamSymmetric(paramSym,BaseAirfoil,calc_gradient):
param = np.r_[paramSym,paramSym[1],-paramSym[2],-paramSym[3],-paramSym[4]]
[delta,dDeltadparam,Shape,Airfoil0] = dDeltadParam(param,BaseAirfoil,calc_gradient)
if calc_gradient:
dParamdParamIn = np.r_[np.eye(5),np.c_[np.zeros((4,1)),np.diag([1,-1,-1,-1])]]
dDeltadparam = dDeltadparam@dParamdParamIn
return [delta,dDeltadparam,Shape,Airfoil0]
else:
return [delta,None,Shape,Airfoil0]
# This function calculates the cross-sectional area of an airfoil (polygon formula)
# We need that later in our optimization function
def AreaPoly(BASE,vd=None,calc_gradient=False):
N = BASE.shape[1]
X0 = BASE[0,:]
Y0 = BASE[1,:]
[Normals,S] = GetNormals(BASE)
dX = np.asarray(Normals[0,:]).flatten()
dY = np.asarray(Normals[1,:]).flatten()
if not vd is None:
XN = X0+vd*(dX)
YN = Y0+vd*(dY)
else:
XN = X0
YN = Y0
Area = 0.5*np.dot((XN[1::]+XN[0:N-1]),(YN[1::]-YN[0:N-1]))
Area_vd = None
if calc_gradient:
Area_vd = np.zeros(N)
Area_vd[1::] = 0.5*(dX[1::])*(YN[1::]-YN[0:N-1]) + 0.5*(dY[1::])*(XN[1::]+XN[0:N-1])
Area_vd[0:N-1] += 0.5*(dX[0:N-1])*(YN[1::]-YN[0:N-1]) - 0.5*(dY[0:N-1])*(XN[1::]+XN[0:N-1])
return (Area,Area_vd,S)
In the following section scipys optimization routine SLSQP(gradient based Newton-type method) is used to optimize the parameters for minimal drag for a symmetric airfoil.
from scipy import optimize
# What is our goal?
# We want to optimize the drag at zero angle of attack (and zero lift) of an airfoil of
# a certain cross-sectional area.
# To that effect we want
# min_param cd
# s.t. A(param) = fac*A0
# where cd is a function of viiflow which may or may not have converged based on the random parameters thrown at it by the optimizer.
# This function is optimized
def evalFun(z,x0,bl,p,setup,BaseAirfoil):
[delta,dDeltadparam,_,_] = dDeltadParamSymmetric(z,BaseAirfoil,True)
_,flag,res,grad,gradients = vf.iter(x0,bl,p,setup,None,None,delta)
cd = bl[0].CD
# If unconverged, let the optimizer know this is not a valid step
if flag<=0 and np.linalg.norm(res)>1e-4:
cd = np.inf
(p,bl,_) = vf.init(BaseAirfoil,setup)
_,_,_,_,gradients = vf.iter(x0,bl,p,setup)
#fig,ax = plt.subplots(1,1)
#ax.plot(res)
fun = cd*1000
grad = 1000*gradients.total.cd_vd[0,:]@dDeltadparam
# We don't really need all these outputs
# But I did need them while debugging the different parts of our function or plotting
return [fun, grad, cd]
def constraintFun(param,A0,BaseAirfoil,calc_gradient):
[delta,dDeltadparam,_,_] = dDeltadParamSymmetric(param,BaseAirfoil,calc_gradient)
# Area
A,A_vd,_ = AreaPoly(BaseAirfoil,delta,calc_gradient)
Adif = (A-A0)/A0
if calc_gradient:
Adif_vd = A_vd/A0
return [Adif,Adif_vd@dDeltadparam]
else:
return [Adif,None]
# This function is called in addition after a successful step
# Here we could also use plot callbacks
def iterfun(z,history,x,bl,p,setup):
history.append(bl[0].CD)
def setup_problem(AreaFac,BaseAirfoil,setup,initial_guess):
# Initial Area
Area0,Area0_vd,S = AreaPoly(BaseAirfoil)
# Our bounds on param for reasonable airfoils
bounds = optimize.Bounds(lb=np.r_[0.001,0.05,-.5,0.01,1*pi/180],
ub=np.r_[.1,.95,-.0001,2,20*pi/180], keep_feasible=False)
# Estimate on Shape
z0 = None
if z0 is None:
# Try to find scaled Naca Parameters
# Is this really necessary? No, but it's fun to see how we could find the parameters to match a shape.
# And the Naca is a good guess, and we don't beed to give initial conditions from previous runs making
# this less reliant on previous convergence or steps taken.
# Trivial guess: Scale what makes sense to scale
z0 = param[0:5].copy()
J = [0,2,3,4]
z0[J] = param[J]*AreaFac
# Newton Loop
if True:
for it in range(20):
# Check the current geometry
[_,_,_,BaseEst] = dDeltadParamSymmetric(z0[0:5],BaseAirfoil,False)
[_,dDeltadparam,_,_] = dDeltadParamSymmetric(z0[0:5],BaseEst,True)
# Scale NACA with Y, find parameter to fit in (one step)
ddelta,_,_,_ = vf.virtual_displacement_from_geometry(BaseEst,np.diag(np.r_[1,AreaFac])@BaseAirfoil)
dz0 = np.linalg.lstsq(dDeltadparam,ddelta,rcond=None)[0]
# The parameters don't change sign, reasonable stepsize limit
lam = min(1.0,np.min(abs(z0)*.2/abs(dz0)))
z0 += lam*dz0
if np.linalg.norm(dz0)<1e-5:
break
[_,_,_,BaseEst] = dDeltadParamSymmetric(z0,BaseAirfoil,False)
BaseEst = vft.repanel(BaseEst,N)
# Setup viiflow and solve NACA0012
(p,bl,x0) = vf.init(BaseEst,setup)
x0,flag,res,grad,_ = vf.iter(x0,bl,p,setup,None,None,None)
constraint = {
'type':'eq',
'fun': lambda z: constraintFun(z,Area0*AreaFac,BaseEst,False)[0],
'jac': lambda z: constraintFun(z,Area0*AreaFac,BaseEst,True)[1]
}
fun = lambda z: evalFun(z,x0,bl,p,setup,BaseEst)[0:2]
if not initial_guess is None:
z0 = initial_guess
return (fun, constraint, z0, bounds,p,bl,setup,x0,BaseEst)
# Set up problem
tol = 1e-5
AreaFac = 1.0
setup = vf.setup(Re=1e6,Ma=0.0,Ncrit=9,Alpha=0,Gradients=True,Silent=True,Tolerance=1e-6)
(fun, constraint, z0, bounds,p,bl,setup,x0,BaseEst) = setup_problem(AreaFac,BaseAirfoil.copy(),setup,initial_guess=None)
method = 'SLSQP'
cd0 = bl[0].CD
history = [cd0]
resultAnalytical = optimize.minimize(fun, z0, method=method, jac= True,bounds=bounds,constraints = constraint ,options={"ftol":tol,"maxiter":200,"disp":True},callback=lambda z: iterfun(z,history,x0,bl,p,setup))
Optimization terminated successfully (Exit mode 0)
Current function value: 4.109520154393415
Iterations: 44
Function evaluations: 105
Gradient evaluations: 44
fig,ax = plt.subplots(1,1)
ax.plot(np.arange(1,len(history)+1),history,'.-',label="SLSQP (Gradient-based)")
#ax.plot(np.arange(1,len(historyHeuristic)+1),historyHeuristic,label="Nelder-Mead")
ax.set_ylabel("Drag coefficient")
ax.set_xlabel("Iteration");
Now we calculate for a range of thicknesses and three Reynolds numbers.
def run_re_updown(varin):
i,BaseAirfoil,tol = varin
REs = [2e5,2e5,5e5,5e5,1e6,1e6]
downs = np.logical_not(np.asarray([True,False,True,False,True,False]))
RE = REs[i]
down = downs[i]
NSteps = 8
results = {}
if down:
AreaFac = np.linspace(1.0,0.5,int(NSteps))
else:
AreaFac = np.linspace(1.0,2,int(NSteps))[1::]
NSteps = len(AreaFac)
DragVector = np.zeros(NSteps)*np.nan
AreaVector = np.zeros(NSteps)*np.nan
XVector = np.zeros((5,NSteps))*np.nan
SuccessVector = np.zeros(NSteps)*np.nan
tol = 1e-5
z0 = None
for k in range(NSteps):
setup = vf.setup(Re=RE,Ma=0.0,Ncrit=9,Alpha=0,Gradients=True,Silent=True,Tolerance=1e-6)
success = True
try:
(fun, constraint, z0, bounds,p,bl,setup,x0,BaseEst) = setup_problem(AreaFac[k],BaseAirfoil,setup,initial_guess=z0)
resultPareto = optimize.minimize(fun, z0, method='SLSQP', jac= True,bounds=bounds,constraints = constraint ,options={"ftol":tol,"maxiter":200})
except:
success = False
print("RE %3.0E\tA/A0 %4.4g\tCD %g\tFlag %i\tfevals %u"%(RE,AreaFac[k], resultPareto.fun/1000,resultPareto.status,resultPareto.nfev))
# If unsuccessful, try again with new guess
if (not resultPareto.success) or (not success) or (k>0 and down and DragVector[k-1]<resultPareto.fun/1000):
(fun, constraint, z0, bounds,p,bl,setup,x0,BaseEst) = setup_problem(AreaFac[k],BaseAirfoil,setup,initial_guess=None)
resultPareto = optimize.minimize(fun, z0, method='SLSQP', jac= True,bounds=bounds,constraints = constraint ,options={"ftol":tol,"maxiter":200})
print("(redo) RE %3.0E\tA/A0 %4.4g\tCD %g\tFlag %i\tfevals %u"%(RE,AreaFac[k], resultPareto.fun/1000,resultPareto.status,resultPareto.nfev))
AreaVector[k] = AreaFac[k]
DragVector[k] = resultPareto.fun/1000
XVector[:,k] = resultPareto.x
SuccessVector[k] = resultPareto.success
if resultPareto.success:
z0 = resultPareto.x
else:
z0 = None
# If the previous one did not converge, retry with this as a starting value
if k>0 and resultPareto.success:
if not SuccessVector[k-1] or (not down and DragVector[k-1]>resultPareto.fun/1000):
(fun, constraint, z01, bounds,p,bl,setup,x0,BaseEst) = setup_problem(AreaFac[k-1],BaseAirfoil,setup,initial_guess=z0)
resultPareto = optimize.minimize(fun, z01, method='SLSQP', jac= True,bounds=bounds,constraints = constraint ,options={"ftol":tol,"maxiter":200})
print("(redo) RE %3.0E\tA/A0 %4.4g\tCD %g\tFlag %i\tfevals %u"%(RE,AreaFac[k-1], resultPareto.fun/1000,resultPareto.status,resultPareto.nfev))
if resultPareto.success and ((DragVector[k-1] > resultPareto.fun/1000) or not SuccessVector[k-1]):
DragVector[k-1] = resultPareto.fun/1000#np.asarray(bl_local[0].CD)
XVector[:,k-1] = resultPareto.x
SuccessVector[k-1] = resultPareto.success
results = {}
results["RE"]=RE
results["Area"]=AreaVector
results["Drag"]=DragVector
results["X"]=XVector
results["Success"] = SuccessVector
return results
from concurrent.futures import ThreadPoolExecutor
with ThreadPoolExecutor(max_workers = 8) as exc:
resultsParallel = exc.map(run_re_updown,[(i,BaseAirfoil.copy(),tol) for i in range(6)])
# Join results
results = list(resultsParallel)
RE 5E+05 A/A0 1.143 CD 0.00597255 Flag 0 fevals 42 RE 5E+05 A/A0 1 CD 0.00564209 Flag 0 fevals 45 RE 2E+05 A/A0 1.143 CD 0.0092399 Flag 0 fevals 62 RE 1E+06 A/A0 1 CD 0.00410976 Flag 0 fevals 57 RE 5E+05 A/A0 0.9286 CD 0.00547719 Flag 0 fevals 35 RE 1E+06 A/A0 1.143 CD 0.00435844 Flag 0 fevals 89 RE 5E+05 A/A0 0.8571 CD 0.00532872 Flag 0 fevals 46 RE 2E+05 A/A0 1 CD 0.00875451 Flag 0 fevals 709 RE 2E+05 A/A0 0.9286 CD 0.00851406 Flag 0 fevals 42 RE 2E+05 A/A0 0.8571 CD 0.00829193 Flag 0 fevals 40 RE 5E+05 A/A0 1.286 CD 0.00629549 Flag 0 fevals 1085 RE 2E+05 A/A0 0.7857 CD 0.00805146 Flag 0 fevals 63 RE 2E+05 A/A0 0.7143 CD 0.00781655 Flag 0 fevals 45 RE 2E+05 A/A0 0.6429 CD 0.00757842 Flag 0 fevals 32 RE 2E+05 A/A0 0.5714 CD 0.00715803 Flag 0 fevals 81 RE 2E+05 A/A0 0.5 CD 0.00686904 Flag 0 fevals 54 RE 2E+05 A/A0 1.286 CD 0.00973548 Flag 9 fevals 2114 RE 5E+05 A/A0 1.429 CD 0.00666086 Flag 0 fevals 522 RE 5E+05 A/A0 0.7857 CD 0.00519281 Flag 0 fevals 1642 RE 1E+06 A/A0 0.9286 CD 0.00399701 Flag 9 fevals 2122 RE 1E+06 A/A0 1.286 CD 0.00448078 Flag 9 fevals 2138 (redo) RE 1E+06 A/A0 1.286 CD 0.00462503 Flag 0 fevals 56 (redo) RE 2E+05 A/A0 1.286 CD 0.00979969 Flag 0 fevals 873 RE 2E+05 A/A0 1.429 CD 0.0102705 Flag 0 fevals 247 RE 5E+05 A/A0 0.7143 CD 0.00502292 Flag 9 fevals 2038 (redo) RE 5E+05 A/A0 0.7143 CD 0.00502124 Flag 0 fevals 395 RE 5E+05 A/A0 0.6429 CD 0.00486489 Flag 0 fevals 30 (redo) RE 1E+06 A/A0 0.9286 CD 0.00399628 Flag 9 fevals 2021 RE 1E+06 A/A0 1.429 CD 0.00462503 Flag 0 fevals 56 (redo) RE 1E+06 A/A0 1.429 CD 0.00487992 Flag 0 fevals 41 RE 5E+05 A/A0 0.5714 CD 0.0047059 Flag 0 fevals 172 RE 5E+05 A/A0 1.571 CD 0.00703039 Flag 9 fevals 2139 RE 5E+05 A/A0 0.5 CD 0.00449611 Flag 0 fevals 74 RE 1E+06 A/A0 1.571 CD 0.00517164 Flag 0 fevals 180 RE 2E+05 A/A0 1.571 CD 0.0106974 Flag 9 fevals 2112 RE 1E+06 A/A0 0.8571 CD 0.00387539 Flag 9 fevals 2004 (redo) RE 2E+05 A/A0 1.571 CD 0.0108307 Flag 9 fevals 2098 RE 1E+06 A/A0 1.714 CD 0.00547163 Flag 0 fevals 1995 RE 1E+06 A/A0 1.857 CD 0.00955117 Flag 0 fevals 122 (redo) RE 5E+05 A/A0 1.571 CD 0.00714789 Flag 9 fevals 2048 RE 5E+05 A/A0 1.714 CD 0.00754648 Flag 0 fevals 58 (redo) RE 5E+05 A/A0 1.571 CD 0.00714281 Flag 0 fevals 187 (redo) RE 1E+06 A/A0 0.8571 CD 0.00387539 Flag 9 fevals 2004 RE 1E+06 A/A0 0.7857 CD 0.00376387 Flag 0 fevals 51 (redo) RE 1E+06 A/A0 0.8571 CD 0.00387585 Flag 0 fevals 99 RE 1E+06 A/A0 0.7143 CD 0.00365835 Flag 0 fevals 44 RE 1E+06 A/A0 0.6429 CD 0.00353611 Flag 0 fevals 42 RE 1E+06 A/A0 0.5714 CD 0.00342664 Flag 0 fevals 36 RE 1E+06 A/A0 0.5 CD 0.00330025 Flag 0 fevals 26 RE 1E+06 A/A0 2 CD 0.00646775 Flag 9 fevals 1504 (redo) RE 1E+06 A/A0 2 CD 0.00621148 Flag 0 fevals 35 RE 2E+05 A/A0 1.714 CD 0.0115851 Flag 0 fevals 1274 RE 5E+05 A/A0 1.857 CD 0.00789299 Flag 9 fevals 2105 (redo) RE 1E+06 A/A0 1.857 CD 0.00601684 Flag 0 fevals 1444 (redo) RE 2E+05 A/A0 1.571 CD 0.0108632 Flag 0 fevals 1228 RE 2E+05 A/A0 1.857 CD 0.0121364 Flag 9 fevals 2056 (redo) RE 5E+05 A/A0 1.857 CD 0.00786551 Flag 9 fevals 2060 RE 5E+05 A/A0 2 CD 0.00835388 Flag 0 fevals 70 (redo) RE 2E+05 A/A0 1.857 CD 0.0121499 Flag 9 fevals 2046 RE 2E+05 A/A0 2 CD 0.012764 Flag 0 fevals 53 (redo) RE 5E+05 A/A0 1.857 CD 0.00793457 Flag 9 fevals 2030 (redo) RE 2E+05 A/A0 1.857 CD 0.0121359 Flag 9 fevals 2100
fig, ax = plt.subplots(1,2)
counter = 0
for k in range(0,len(results),2):
Areas = np.r_[np.asarray(results[k+1]["Area"])[::-1],np.asarray(results[k]["Area"])]
Drags = np.r_[np.asarray(results[k+1]["Drag"])[::-1],np.asarray(results[k]["Drag"])]
X = np.c_[np.asarray(results[k+1]["X"])[:,::-1],np.asarray(results[k]["X"])]
RE = results[k]["RE"]
for j in [0,int(len(Areas)/2),len(Areas)-1]:
off = -.2+.4*(j/(len(Areas)-1))
[_,_,_,Shape] = dDeltadParamSymmetric(X[:,j],BaseAirfoil,False)
ax[1].plot(Shape[0,:],Shape[1,:]+off,color = "C%u"%counter)
ax[0].plot(Areas,Drags,'.-',label="Optimum RE %1.1e"%RE)
# Add non-converged hint
#for j in range(len(results[k]["Area"])):
# if not results[k]["Success"][j]:
# ax[0].plot(results[k]["Area"][j],results[k]["Drag"][j],'x')
ax[1].set_prop_cycle(None)
counter+=1
ax[0].set_prop_cycle(None)
# Add simple NACA solution
for RE in [2e5,5e5,1e6]:
setup.Re = RE
(p,bl,x) = vf.init(BaseAirfoil,setup)
vf.iter(x,bl,p,setup)
ax[0].plot(1.0,bl[0].CD,'v',label="NACA0012 RE %1.1e"%RE)
ax[1].text(.3, .2, '200%', fontsize = 12,ha='center', va='center')
ax[1].text(.3, 0, '100%', fontsize = 12,ha='center', va='center')
ax[1].text(.3, -.2, '50%', fontsize = 12,ha='center', va='center')
#ax[1].plot(BaseAirfoil[0,:],BaseAirfoil[1,:],'--k',label='NACA0012')
#ax[1].set_aspect('equal', adjustable='box')
ax[0].set_ylim(bottom=0.0)
ax[0].set_xlabel('Area/NACA0012')
ax[0].set_ylabel('CD')
ax[1].set_xlabel('x/c')
ax[0].legend()
<matplotlib.legend.Legend at 0x1118960f0>
The higher the Reynolds Number, the smaller the leading-edge radius and the more area is put aft - or the less "raindroppy" the airfoil becomes. The improvement of drag from the starting airfoil is smallest at 5e5, where the NACA0012 seems to be close to the optimum shape for minimum drag at these settings.
[1] Drela, Mark. Pros and cons of airfoil optimization. Frontiers of computational fluid dynamics 1998 (1998): 363-381.
[2] Sobieczky, Helmut. Parametric airfoils and wings. Recent development of aerodynamic design methodologies. Vieweg+ Teubner Verlag, 1999. 71-87.
[3] Yu, J. C., and R. Wulandari. Airfoil aerodynamics optimization under uncertain operating conditions. Journal of Physics: Conference Series. Vol. 1446. No. 1. IOP Publishing, 2020.
[4] Della Vecchia, Pierluigi, Elia Daniele, and Egidio DʼAmato. An airfoil shape optimization technique coupling PARSEC parameterization and evolutionary algorithm. Aerospace Science and Technology 32.1 (2014): 103-110.
[5] Drela, Mark. XFOIL: An analysis and design system for low Reynolds number airfoils. Low Reynolds number aerodynamics. Springer, Berlin, Heidelberg, 1989. 1-12.