Generally speaking, for low-speed flow around an aerofoil, the largest gradients in the flow-properties coincide with the greatest curvature of the aerofoil, and thus are the region that requires the greatest concentration of cells (or nodes) in the computational domain.
Many of the simplified treatments for spacing points along an aerofoil use a logarithmic spacing biased towards the leading edge. This of course can create issues if spaced in the chord-wise direction (x/c). The surface-normal at the leading edge of a symmetric aerofoil is perpendicular to the stream-wise direction, and thus can become under-resolved, and even faceted, with a simple linear or logarithmic chord-wise spacing.
One of the more significant issues with linear or logarithmic spacing is that it simply does not discriminate on the basis of curvature, thus the growth factor becomes a qualitative parameter controlled by the aerodynamicist. Any other curvature, for example the upper-surface curvature caused by the s-shaped chord of a DHTMU aerofoil, is neglected.
Some computational treatments, though able to nicely capture curvature, are a little too complex for this simplified tool intended for 2D aerofoil analyses. I believe that the Ramer-Douglas-Peucker (RDP) algorithm is a simple yet powerful tool ideally suited to be applied to spacing points/nodes along a computational grid for for 2D aerofoil analyses.
RDP is a a decimation (down-sampling) algorithm, which although cannot increase resolution (create points) in areas of high-curvature, it can remove remove points in areas where curvature is low. This approach is quite simple and elegant. The algorithm works by creating a line between the start and end points, then recursively divides this line whilst each time finding the furthest point, marking these points to be kept if they are greater than a user-specified distance (epsilon).
As this requires a set of points, rather than a spline or other analytical curve, there are upsides and downsides. In the case of of analytic aerofoils, such as the DHTMU, or the NACA 4 that I originally scripted into GEVgoil, these need to be discretised before being down-sampled. However, for aerofoils that are already being supplied as a set of points, for example using the CSV method that is now also in GEVfoil, RDP can be applied directly.
To simplify matters further, RDP is available as a Python module via PIP: https://pypi.org/project/rdp/ Using this excellent ready-made library, I have constructed perhaps the simplest Python class, implementing this RDP module into GEVfoil:
from rdp import rdp
import numpy as np
"""A class to decimate aerofoil points"""
def __init__(self, x, y, e=0.0001):
xy = np.vstack((x, y)).T
self.reducedXY = rdp(xy, epsilon=e)
self.reducedX = self.reducedXY[:,0]
self.reducedY = self.reducedXY[:,1]
And of course, the proof is in the pudding. Some examples of the RDP module in action in GEVfoil with varying epsilon:
Previously I had made a simple function that will generate a number of points/coordinates for any specified NACA 4 series aerofoil. The limitation with this approach was that it did a fairly poor job at capturing the leading-edge. The leading edge is of critical importance to the aerodynamic performance and characteristics of the aerofoil. After some consideration, I have decided not to modify this technique. I am designing the GEVfoil tools to be as modular, generic and extendable as possible. This isn’t to say I am neglecting the importance of the leading edge, or that the method is flawed.
Initially I had considered returning the aerofoil as a mathematical function, rather than a set of points. This concept may well be revisited in the future. However, most of the other input aerofoils will arrive as a set of coordinates of varying fidelity, from sources such as a the UIUC Airfoil Coordinates Database or from Comma Separated Variable (CSV) files. Therefore, the classes which will receive the aerofoil coordinate objects will deal with interpolating and appropriately converting coordinates into functions anyway. Consequently, as long as the linear-spacing is arbitrarily high enough a simple piece-wise cubic interpolation will suffice.
Note that the CSV reader itself is an abstraction of the AerofoilReader class, which I plan to later extend to read the DAT file from UIUC,
To that end, the CSV reader currently looks as follows:
""" An abstract class to be inherited by classes for specific
aerofoil file types """
def __init__(self, filename):
self.coord = 
self.title =  # Column titles
self.label =  # Aerofoil Label
foil = open(filename, 'r')
self.lines = foil.readlines()
""" Read in CSV. Assumes a header is present."""
def __init__(self, filepath, xcol, ycol, zcol=False):
self.label = os.path.basename(filepath).split('.')
inpcols = 3 # Two or three columns
if not zcol:
inpcols = 2
for i, line in enumerate(self.lines):
parts = line.rstrip().replace('"', '').split(',')
additional = 
if i > 0:
# Not the header
if len(parts) > inpcols:
# More than x, y & z - store other columns
for j in range(len(parts)):
if j != xcol and j != ycol and j != zcol:
xyz = [parts[xcol], parts[ycol], parts[zcol]]
self.coord.append(xyz + additional)
# It's the header
self.title = parts
The intention behind GEVfoil is to establish an open-source set of Python classes to handle aerofoils (American English: airfoils). The intention is that these can be used for anything from generating simple point clouds, analytical representations all the way to STL files for CFD meshing and OpenFOAM BlockMesh output.
The first step is starting with something basic, in this case a simple class to generate NACA 4 series aerofoils:
import numpy as np
import matplotlib.pyplot as plt
"""A class for generating NACA 4 series aerofoils
TO-DO - linspace not good at LE capture """
def __init__(self, foil, points=200, chord=1):
m = float(foil)/100 # max camber
p = float(foil)/10 # chordwise position of max camber
t = float(foil[2:])/100 # thickness
x = np.linspace(0, chord, points)
self.x = x
yt = self.thickness(x, t)
yc = 
for coord in x:
yc.append(self.camber(m, p, coord))
# No camber
y1 = yc + yt
y2 = yc - yt
self.y = y1
plt.plot(self.x, y1, '.')
plt.plot(self.x, y2, '.')
def thickness(self, x, t):
# Y thickness at given x point
return t / 0.2 * \
def camber(self, m, p, x):
# Return the camber of the aerofoil
if x <= p:
return m/p**2 * (2*p*x-x**2)
if __name__ == "__main__":
test = NACA4('2412')
For now, Matplotlib has been built into the class for visual reference. For example, 2412 generates:
And 0010 creates:
At this stage there is a clear issue around the leading edge, the function was interrogated with a linearly spaced x/chord-wise interval, this must be improved for the future to cluster points around the greatest curvature.