"""
.. airfoil
-------
Contains a class for creating, modifying and exporting airfoils.
"""
import numpy as np
from pyspline import Curve
from scipy.optimize import brentq, newton, minimize
from . import sampling
from .utils.geom_ops import _translateCoords, _rotateCoords, _scaleCoords, _getClosestY
from .utils.io import _writeDat, _writePlot3D, _writeFFD, Error
EPS = np.finfo(np.float64).eps
ZEROS_2 = np.zeros(2)
[docs]class Airfoil:
"""
A class for manipulating airfoil geometry.
Create an instance of an airfoil. There are two ways of instantiating
this object: by passing in a set of points, or by reading in a coordinate
file. The points must satisfy the following requirements:
- Ordered such that they form a continuous airfoil surface
- First and last points correspond to trailing edge
- If the two points coincide, the airfoil will be considered sharp/round
- Otherwise, a blunt (open) TE will be assumed
It is not necessary for the points to be in a counter-clockwise ordering. If
they are not ordered counter-clockwise, the order will be reversed so that
all functions can be written to expect the spline to begin at the upper
surface of the trailing edge and end at the lower surface of the trailing
edge.
See documentation in :class:`pySpline <pyspline:pyspline.pyCurve.Curve>` for information on the spline representation
Parameters
----------
coords : ndarray[N,3]
Full array of airfoil coordinates
spline_order : {4, 2, 3}
Order of the spline. :math:`n` order implies :math:`C^{n-2}` continuity
normalize : bool
True to normalize the chord of the airfoil, set to zero angle of attack, and move the leading edge to the origin
nCtl : int
If this is set to an integer the underlying airfoil spline will be created using least mean squares (LMS) instead of interpolation. The number of control points used in the LMS will be equal to nCtl.
"""
def __init__(self, coords, spline_order=4, normalize=False, nCtl=None):
self.spline_order = spline_order
self.sampled_pts = None
self.closedCurve = None
self.camber = None
self.british_thickness = None
self.american_thickness = None
self.nCtl = nCtl
# Initialize geometric information
self.recompute(coords)
if normalize:
self.normalizeAirfoil()
[docs] def recompute(self, coords):
"""
Recomputes the underlying spline and relevant parameters from the given set of coordinates.
Parameters
----------
coords : Ndarray [N,2]
The coordinate pairs to compute the airfoil spline from
"""
if self.nCtl:
self.spline = Curve(X=coords, k=self.spline_order, nCtl=self.nCtl)
else:
self.spline = Curve(X=coords, k=self.spline_order)
self.reorder()
self.TE = self.getTE()
self.LE, self.s_LE = self.getLE()
self.chord = self.getChord()
self.twist = self.getTwist()
self.closedCurve = (coords[0, :] == coords[-1, :]).all()
self.sampled_pts = None
camber_pts = self.getCDistribution(coords.size)
self.camber = Curve(X=camber_pts, k=3)
self.british_thickness = Curve(X=self.getThickness(coords.size, "british"), k=3)
self.american_thickness = Curve(X=self.getThickness(coords.size, "american"), k=3)
[docs] def reorder(self):
"""
This function orients the points counterclockwise and sets the start point to the TE
"""
# Check to make sure spline ends at TE (For now assume this is True)
# Make sure oriented in counter-clockwise direction.
coords = self.spline.X
N = coords.shape[0]
orientation = 0
for i in range(1, N - 1):
v = coords[i + 1] - coords[i]
r = coords[i + 1] - coords[i - 1]
# skip duplicate points
if np.linalg.norm(r) < EPS:
continue
s = (coords[i, 0] * r[0] + coords[i, 1] * r[1]) / np.linalg.norm(r)
n = coords[i] - r * s
if np.linalg.norm(n) > EPS:
n = n / np.linalg.norm(n)
orientation += n[0] * v[1] - n[1] * v[0]
if orientation < 0:
# Flipping orientation to counter-clockwise
self.recompute(self.spline.X[::-1, :])
## Geometry Information
[docs] def getCamber(self):
"""
Calculates the camber spline defined by the airfoil
Returns
-------
camber : pySpline curve object
The spline that defines the camberline from s = 0 at the leading edge to s = 1 at the trailing edge.
"""
return self.camber
[docs] def getTE(self):
"""
Calculates the trailing edge point on the spline
Returns
-------
TE : 2darry [2]
The coordinate of the trailing edge of the airfoil
"""
TE = (self.spline.getValue(0) + self.spline.getValue(1)) / 2
return TE
[docs] def getLE(self):
r"""
Calculates the leading edge point on the spline, which is defined as the point furthest away from the TE. The spline is assumed to start at the TE. The routine uses a root-finding algorithm to compute the LE.
Returns
-------
LE : Ndarray [2]
the coordinate of the leading edge
s_LE : float
the parametric position of the leading edge
Notes
-----
Let the TE be at point :math:`x_0, y_0`, then the Euclidean distance between the TE and any point on the airfoil spline is :math:`\ell(s) = \sqrt{\Delta x^2 + \Delta y^2}`, where :math:`\Delta x = x(s)-x_0` and :math:`\Delta y = y(s)-y_0`. We know near the LE, this quantity is concave. Therefore, to find its maximum, we differentiate and use a root-finding algorithm on its derivative.
:math:`\\frac{\mathrm{d}\ell}{\mathrm{d}s} = \\frac{\Delta x\\frac{\mathrm{d}x}{\mathrm{d}s} + \Delta y\\frac{\mathrm{d}y}{\mathrm{d}s}}{\ell}`
The function ``dellds`` computes the quantity :math:`\Delta x\\frac{\mathrm{d}x}{\mathrm{d}s} + \Delta y\\frac{\mathrm{d}y}{\mathrm{d}s}` which is then used by ``brentq`` to find its root, with an initial bracket at :math:`[0.3, 0.7]`.
"""
def dellds(s, spline, TE):
pt = spline.getValue(s)
deriv = spline.getDerivative(s)
dx = pt[0] - TE[0]
dy = pt[1] - TE[1]
return dx * deriv[0] + dy * deriv[1]
s_LE = brentq(dellds, 0.3, 0.7, args=(self.spline, self.TE))
LE = self.spline.getValue(s_LE)
return LE, s_LE
[docs] def getTwist(self):
"""
Calculates the twist of the airfoil using the leading and trailing edge points
Returns
-------
twist : float
The twist in degrees of the airfoil
"""
chord_vec = self.TE - self.LE
twist = np.arctan2(chord_vec[1], chord_vec[0]) * 180 / np.pi
# twist = np.arccos(normalized_chord.dot(np.array([1., 0.]))) * np.sign(normalized_chord[1])
return twist
[docs] def getChord(self):
"""
Calculates the chord of the airfoil as the distance between the leading and trailing edges
Returns
-------
chord : float
The chord length
"""
chord = np.linalg.norm(self.TE - self.LE)
return chord
[docs] def getSplinePts(self):
"""
alias for returning the points that make the airfoil spline
Returns
-------
X : Ndarry [N, 2]
the coordinates that define the airfoil spline
"""
return self.spline.X
[docs] def findPt(self, position, axis=0, s_0=0):
"""
finds that point at the intersection of the plane defined by the axis and the postion
and the airfoil curve
Parameters
----------
position : float
the position of the plane on the given axis
axis : int
the axis the plane will intersect 0 for x and 1 for y
s_0 : float
an initial guess for the parameteric position of the solution
Returns
-------
X : Ndarray [2]
The coordinate at the intersection
s_x : float
the parametric location of the intersection
"""
def err(s):
return self.spline(s)[axis] - position
def err_deriv(s):
return self.spline.getDerivative(s)[axis]
s_x = newton(err, s_0, fprime=err_deriv)
return self.spline.getValue(s_x), s_x
[docs] def getTEThickness(self):
"""
gets the trailing edge thickness for the airfoil
Returns
-------
TE_thickness : float
the trailing edge thickness
"""
top = self.spline.getValue(0)
bottom = self.spline.getValue(1)
TE_thickness = np.sqrt((top[0] - bottom[0]) ** 2 + (top[1] - bottom[1]) ** 2)
return TE_thickness
[docs] def getLERadius(self):
"""
Computes the leading edge radius of the airfoil. Note that this is heavily dependent on the initialization points, as well as the spline order/smoothing.
Returns
-------
LE_rad : float
The leading edge radius
"""
# if self.s_LE is None:
# self.getLE()
first = self.spline.getDerivative(self.s_LE)
second = self.spline.getSecondDerivative(self.s_LE)
LE_rad = np.linalg.norm(first) ** 3 / np.linalg.norm(first[0] * second[1] - first[1] * second[0])
return LE_rad
[docs] def getCDistribution(self, nPts):
"""
Return the coordinates of the camber points
Parameters
----------
nPts : int
The number of points to sample
Returns
-------
camber_pts : Ndarray [nPts, 2]
the locations of the camber points of the airfoil starting with the leading edge and ending with the trailing edge
"""
top_surf, bottom_surf = self.splitAirfoil()
# Compute the chord
chord_pts = np.vstack([self.LE, self.TE])
chord = Curve(X=chord_pts, k=2)
# Sampling along airfoil for camber points
lin_sampling = np.linspace(0, 1, nPts - 1, endpoint=False)[1:]
chord_pts = chord.getValue(lin_sampling)
camber_pts = np.zeros((nPts - 2, 2))
# At each point we are looking for the camber
for j in range(chord_pts.shape[0]):
# Get the direction normal to the chord line
direction = np.array(
[np.cos(np.pi / 2 + np.deg2rad(self.twist)), np.sin(np.pi / 2 + np.deg2rad(self.twist))]
)
direction = direction / np.linalg.norm(direction)
# Draw a ray through the airfoil in the given direction
top = chord_pts[j, :] + 1 * self.chord * direction
bottom = chord_pts[j, :] - 1 * self.chord * direction
temp = np.vstack((top, bottom))
normal = Curve(X=temp, k=2)
# Determine the intersection of this ray with both the upper and lower surfaces
s_top, _, _ = top_surf.projectCurve(normal, nIter=5000, eps=EPS)
s_bottom, _, _ = bottom_surf.projectCurve(normal, nIter=5000, eps=EPS)
intersect_top = top_surf.getValue(s_top)
intersect_bottom = bottom_surf.getValue(s_bottom)
# Compute the camber
camber_pts[j, :] = (intersect_top + intersect_bottom) / 2
# Add TE and LE to the camber points.
camber_pts = np.vstack((self.LE, camber_pts, self.TE))
return camber_pts
[docs] def getThickness(self, nPts, tType):
"""
Computes the thicknesses at each x stations spaced linearly along the airfoil
Parameters
----------
nPts : int
number of points to sample including the edge
tType : str
either "american" or "british"
Returns
-------
thickness_pts : Ndarray [nPts, 2]
The thickness at each x station
"""
if tType not in ["american", "british"]:
raise Error("Do not recognize thickness type!")
top_surf, bottom_surf = self.splitAirfoil()
# The parametric spline values along the camber line to find thickness points
s = np.linspace(0, 1, nPts - 1, endpoint=False)[1:]
thickness_pts = np.zeros((nPts - 2, 2))
# Find thickness at each point
for j in range(len(s)):
# If british we project a ray normal to chordline
if tType == "british":
direction = np.array(
[np.cos(np.pi / 2 - np.deg2rad(self.twist)), np.sin(np.pi / 2 - np.deg2rad(self.twist))]
)
# If american we project a ray normal to camberline
else:
dx = self.camber.getDerivative(s[j])
direction = np.array([-dx[1], dx[0]])
direction = direction / np.linalg.norm(direction)
# create a ray through the upper and lower surfaces from given direction
top = self.camber.getValue(s[j]) + 10 * self.chord * direction
bottom = self.camber.getValue(s[j]) - 10 * self.chord * direction
normal = Curve(X=np.vstack([top, bottom]), k=2)
# approximate location of the intersection as a percentage of the chord
s_guess = (normal.getValue(0.5)[0] - self.LE[0]) / self.chord
# top surf goes from 0 at TE to 1 at LE, so parameter needs to be reversed
top_guess = 1 - s_guess
bottom_guess = s_guess
# Keep guesses within the bounds of the spline
if s_guess > 1:
top_guess = 0
bottom_guess = 1
elif s_guess < 0:
top_guess = 1
bottom_guess = 0
# Find upper and lower intersections
s_top, _, _ = top_surf.projectCurve(normal, nIter=100, eps=EPS, s=top_guess, t=0.5)
s_bottom, _, _ = bottom_surf.projectCurve(normal, nIter=100, eps=EPS, s=bottom_guess, t=0.5)
# Compute the thickness
thickness_pts[j, 0] = self.camber.getValue(s[j])[0]
if tType == "british":
thickness_pts[j, 1] = top_surf.getValue(s_top)[1] - bottom_surf.getValue(s_bottom)[1]
else:
x_top = top_surf.getValue(s_top)
x_bottom = bottom_surf.getValue(s_bottom)
thickness_pts[j, 1] = np.linalg.norm(x_top - x_bottom)
# Add the trailing and leading edge points when we return
return np.vstack([[self.LE[0], 0], thickness_pts, [self.TE[0], self.getTEThickness()]])
[docs] def getTEAngle(self):
"""
Computes the trailing edge angle of the airfoil. We assume here that the spline goes from top to bottom, and that s=0 and s=1 corresponds to the
top and bottom trailing edge points. Whether or not the airfoil is closed is irrelevant.
Returns
-------
TE_angle : float
The angle of the trailing edge in degrees
"""
top = self.spline.getDerivative(0)
top = top / np.linalg.norm(top)
bottom = self.spline.getDerivative(1)
bottom = bottom / np.linalg.norm(bottom)
# print(np.dot(top,bottom))
TE_angle = np.pi - np.arccos(np.dot(top, bottom))
return np.rad2deg(TE_angle)
[docs] def getMaxThickness(self, tType):
"""
Parameters
----------
tType : str
Can be one of 'british' or 'american'
Returns
-------
x_loc : float
The x station containing the maximum thickness
max_thickness : float
the maximum thickness of the airfoil
"""
if tType not in ["american", "british"]:
raise Error("Do not recognize thickness type!")
def american_f(s):
return -self.american_thickness.getValue(s)[1]
def american_df(s):
return -self.american_thickness.getDerivative(s)[1]
def british_f(s):
return -self.british_thickness.getValue(s)[1]
def british_df(s):
return -self.british_thickness.getDerivative(s)[1]
if tType == "american":
opt = minimize(american_f, 0.5, method="SLSQP", jac=american_df, bounds=[(0, 1)])
if not opt.success:
raise Error("Could not determine the maximum thickness.")
opt_point = self.american_thickness.getValue(opt.x)
else:
opt = minimize(british_f, 0.5, method="SLSQP", jac=british_df, bounds=[(0, 1)])
if not opt.success:
raise Error("Could not determine the maximum thickness.")
opt_point = self.british_thickness.getValue(opt.x)
return opt_point[0], opt_point[1]
def _findChordProj(self, coord):
"""
Finds the point on the chordline that defines a line from `coord` to the chordline that is perpendicular to the chordline
Parameters
----------
coord : 2darray
The point of interest we wish to project onto the chordline.
Returns
-------
point : 2darray
The coordinate that is the perpendicular projection of `coord` onto the chordline
"""
# vector defines the chord
chord = self.LE - self.TE
# Parametric position of point on chordline
s = (-chord[0] * (self.TE[0] - coord[0]) - chord[1] * (self.TE[1] - coord[1])) / (chord[0] ** 2 + chord[1] ** 2)
return self.TE + s * chord
def _MaxCamberOptimize(self, maximum):
"""
Used to compute the most negative and most positive cambers of an airfoil
Parameters
----------
maximum : bool
If true find most positive, if false find most negative
Returns
-------
x_loc : float
the x location of the maximum camber
max_camber : float
the maximum camber
"""
def f(s, factor):
# Find the perpindicular project onto the chord line
pointInterest = self.camber.getValue(s)
chordProj = self._findChordProj(pointInterest)
# Determine if distance is +/- with cross product
chord = self.LE - self.TE
chord /= np.linalg.norm(chord)
direction = pointInterest - chordProj
direction /= np.linalg.norm(direction)
cross = direction[0] * chord[1] - direction[1] * chord[0]
return cross * factor * np.linalg.norm(pointInterest - chordProj)
if maximum:
factor = -1
else:
factor = 1
opt = minimize(lambda s: f(s, factor), 0.5, method="SLSQP", bounds=[(0, 1)])
if not opt.success:
if maximum:
raise Error("Could not determine maximum camber.")
raise Error("Could not determine minimum camber.")
opt_point = self.camber.getValue(opt.x)
opt_int = self._findChordProj(opt_point)
# convert to airfoil coordinates
x = np.linalg.norm(opt_int - self.LE) / np.linalg.norm(self.LE - self.TE)
c = factor * f(opt.x, factor) / np.linalg.norm(self.LE - self.TE)
return x, c
[docs] def getMaxCamber(self):
"""
Returns
-------
x_loc : float
the x location of the maximum camber
max_camber : float
the maximum camber of the airfoil
"""
return self._MaxCamberOptimize(True)
[docs] def getMinCamber(self):
"""
Returns
-------
x_loc : flaot
the x location of the maximum ngative camber
min_camber : float
the maximum negative camber of the airfoil
"""
return self._MaxCamberOptimize(False)
[docs] def isReflex(self):
"""
Determines if an airfoil is reflex by checking if the derivative of the camber line at the trailing edge is positive.
Returns
-------
reflexive : bool
True if reflexive
"""
if self.camber is None:
self.getCamber()
return self.camber.getDerivative(1)[1] > 0
[docs] def isSymmetric(self, tol=1e-6):
"""
Checks if an airfoil is symmetric
Parameters
----------
tol : float
tolerance for camber line to still be consdiered symmetrical
Returns
-------
symmetric : bool
True if the airfoil is symmetric within the given tolerance
"""
if abs(self.getMinCamber()[1]) < tol and abs(self.getMaxCamber()[1]) < tol:
return True
return False
# ==============================================================================
# Geometry Modification
# ==============================================================================
[docs] def rotate(self, angle, origin=ZEROS_2):
"""
rotates the airfoil about the specified origin
Parameters
----------
angle : float
the angle to rotate the airfoil in degrees
origin : Ndarray [2]
the point about which to rotate the airfoil
"""
new_coords = _rotateCoords(self.spline.X, np.deg2rad(angle), origin)
self.recompute(new_coords)
[docs] def derotate(self, origin=ZEROS_2):
"""
derotates the airfoil about the origin by the twist
Parameters
----------
origin : Ndarray [2]
the location about which to preform the rotation
"""
self.rotate(-1.0 * self.twist, origin=origin)
[docs] def scale(self, factor, origin=ZEROS_2):
"""
Scale the airfoil by factor about the origin
Parameters
----------
factor : float
the scaling factor
origin : Ndarray [2]
the coordinate about which to preform the scaling
"""
new_coords = _scaleCoords(self.spline.X, factor, origin)
self.recompute(new_coords)
[docs] def normalizeChord(self, origin=ZEROS_2):
"""
Set the chord to 1 by scaling the airfoil about the given origin
Parameters
----------
origin : Ndarray [2]
the point about which to scale the airfoil
"""
if self.chord != 1:
self.scale(1.0 / self.chord, origin=origin)
[docs] def translate(self, delta):
"""
Translate the airfoil by the vector delta
Parameters
----------
delta : Ndarray [2]
the vector that defines the translation of the airfoil
"""
coords = _translateCoords(self.spline.X, delta)
self.recompute(coords)
[docs] def center(self):
"""
Move the airfoil so that the leading edge is at the origin
"""
if not np.all(self.LE == ZEROS_2):
self.translate(-1.0 * self.LE)
[docs] def splitAirfoil(self):
"""
Splits the airfoil into upper and lower surfaces
Returns
-------
top : pySpline curve object
A spline that defines the upper surface
bottom : pySpline curve object
A spline that defines the lower surface
"""
# if self.s_LE is None:
# self.getLE()
top, bottom = self.spline.splitCurve(self.s_LE)
return top, bottom
[docs] def normalizeAirfoil(self, derotate=True, normalize=True, center=True):
"""
Sets the twist to zero, the chord to one, and the leading edge location to the origin
Parameters
----------
derotate : bool
True to set twist to zero
normalize : bool
True to set the chord length to one
center : bool
True to put the leading edge at the origin
"""
if derotate or normalize or center:
# Order of operation here is important, even though all three operations are linear, because
# we rotate about the origin for simplicity.
if center:
self.center()
if derotate:
self.derotate()
if normalize:
self.normalizeChord()
[docs] def makeBluntTE(self, xCut=0.98):
"""
This cuts the upper and lower surfaces to creates a blunt trailing edge perpendicular to the chord line.
Parameters
----------
xCut : float
the location to cut the blunt TE **as a percentage of the chord**
"""
# Find global coordinates of cut point
xCut = self.LE + xCut * (self.TE - self.LE)
# The direction normal to the chordline
direction = np.array([np.cos(np.pi / 2 + np.deg2rad(self.twist)), np.sin(np.pi / 2 + np.deg2rad(self.twist))])
direction = direction / np.linalg.norm(direction)
# ray to intersect upper and lower surfaces
ray = [xCut - 2 * direction * self.getChord(), xCut + 2 * direction * self.getChord()]
top_surf, bottom_surf = self.splitAirfoil()
normal = Curve(X=ray, k=2)
# Get intersections
s_top, _, _ = top_surf.projectCurve(normal, nIter=5000, eps=EPS)
s_bottom, _, _ = bottom_surf.projectCurve(normal, nIter=5000, eps=EPS)
# Get all the coordinates that will not be cut off
coords = [top_surf.getValue(s_top)]
chord = self.LE - self.TE
for x in self.getSplinePts():
# dot product test checks for positive projection onto chord
current_direction = x - xCut
if chord[0] * current_direction[0] + chord[1] * current_direction[1] > 0:
coords.append(np.array(x))
coords.append(np.array(bottom_surf.getValue(s_bottom)))
self.recompute(np.array(coords))
[docs] def sharpenTE(self, xCut=0.98):
"""
this method creates a sharp trailing edge **from a blunt one** by projecting straight lines from the upper and lower surfacs of a blunt trailing edge.
Parameters
----------
xCut : float
x location **as a percentage of chord** to cut off the current trailing edge if it is not already blunt.
"""
if xCut >= 1.0 or xCut <= 0:
raise Error("xCut must be between 0 and 1.")
if not self.closedCurve:
self.makeBluntTE(xCut)
# Value of blunt TE point on upper surface
val_u = self.spline.getValue(0)
# derivative of blunt TE point of upper surface wrt parametric parameter
ds_u = self.spline.getDerivative(0)
# slope of blunt TE point of upper surface
dx_u = ds_u[1] / ds_u[0]
# Value of blunt TE point on lower surface
val_l = self.spline.getValue(1)
# derivative of blunt TE point of lower surface wrt parameteric parameter
ds_l = self.spline.getDerivative(1)
# slope of blunt TE point of lower surface
dx_l = ds_l[1] / ds_l[0]
# make sure that the slope of the lower surface is greater than the upper, ensuring the points will intersect
if dx_u == dx_l:
raise Error("Slopes at blunt TE are parallel, no intersection point for a sharp TE.")
if dx_u > dx_l:
raise Error("Slopes at blunt TE indicate an intersection towards the LE of the airfoil.")
# calculate the x location of the intersection
x = (val_l[1] - val_u[1] - val_l[0] * dx_l + val_u[0] * dx_u) / (dx_u - dx_l)
# calculate the y location of the intersection
y = val_l[1] + dx_l * (x - val_l[0])
# add intersection points and then recompute the airfoil
coords = np.vstack(([x, y], self.spline.X, [x, y]))
self.recompute(coords)
[docs] def roundTE(self, xCut=0.98, k=4, nPts=20, dist=0.4):
"""
this method creates a smooth round trailing edge **from a blunt one** using a spline. If the trailing edge is not already blunt xCut specifies the location of the cut
Parameters
----------
xCut : float
x location of the cut **as a percentage of the chord**. Will not do anything if the TE is already blunt.
k: int (3 or 4)
order of the spline used to make the rounded trailing edge of the airfoil.
nPts : int
Number of trailing edge points to add to the airfoil spline
dist : float
Arbitrary factor that specifies how long to make the added TE. Larger dist corresponds to longer addition to the end
"""
if xCut >= 1.0 or xCut <= 0:
raise Error("xCut must be between 0 and 1.")
if not self.closedCurve:
self.makeBluntTE(xCut)
# unit length for making rounded TE
dx = self.getTEThickness() * dist
# create the knot vector for the spline
t = [0] * k + [0.5] + [1] * k
# create the vector of control points for the spline
coeff = np.zeros((k + 1, 2))
for ii in [0, -1]:
coeff[ii] = self.spline.getValue(np.abs(ii))
dX_ds = self.spline.getDerivative(np.abs(ii))
dy_dx = dX_ds[1] / dX_ds[0]
# the indexing here is a bit confusing.ii = 0 -> coeff[1] and ii = -1 -> coef[-2]
coeff[3 * ii + 1] = np.array([coeff[ii, 0] + dx * 0.5, coeff[ii, 1] + dy_dx * dx * 0.5])
if k == 4:
chord = self.TE - self.LE
chord /= np.linalg.norm(chord)
coeff[2] = np.array([self.TE[0] + chord[0] * dx, self.TE[1] + chord[1] * dx])
## make the TE curve
te_curve = Curve(t=t, k=k, coef=coeff)
# ----- combine the TE curve with the spline curve -----
upper_curve, lower_curve = te_curve.splitCurve(0.5)
upper_pts = upper_curve.getValue(np.linspace(1, 0, nPts // 2))
lower_pts = lower_curve.getValue(np.linspace(1, 0, nPts // 2))
coords = np.vstack((upper_pts[:-1], self.spline.X, lower_pts[1:]))
# ---- recompute with new TE ---
self.recompute(coords)
[docs] def removeTE(self, tol=0.3, xtol=0.9):
"""
Removes points from the trailing edge of an airfoil, and recomputes the underlying spline.
Parameters
----------
tol : float
A point is part of the trailing edge if the magnitude of the dot product of the normalized vector of `coord[i+1]-coord[i]` and the normalized vector from the trailing edge to the leading edge is less than this tolerance.
This means that decreasing the tolerance will require the orientation of an element to approach being perpendicular to the chord to be consider part of the trailing edge.
xtol : float
Only checks for trailing edge points if the coodinate is past this fraction of the chord.
Returns
-------
TE_points : Ndarray [N,2]
The points that were flagged as trailing edge points and removed from the airfoil coordinates.
"""
coords = self.getSplinePts()
chord_vec = self.TE - self.LE
unit_chord_vec = chord_vec / np.linalg.norm(chord_vec)
airfoil_mask = []
TE_mask = []
for ii in range(coords.shape[0] - 1): # loop over each element
if coords[ii, 0] >= (self.LE + chord_vec * xtol)[0]:
delta = coords[ii + 1] - coords[ii]
unit_delta = delta / np.linalg.norm(delta)
if np.abs(np.dot(unit_chord_vec, unit_delta)) < tol:
TE_mask += [ii, ii + 1]
else:
airfoil_mask += [ii, ii + 1]
else:
airfoil_mask += [ii, ii + 1]
# list(set()) removes the duplicate pts
self.recompute(coords[list(set(airfoil_mask))])
return coords[list(set(TE_mask))]
## Sampling
[docs] def getSampledPts(self, nPts, spacingFunc=sampling.polynomial, func_args=None, nTEPts=0, TE_knot=False):
"""
This function defines the point sampling along the airfoil surface.
Parameters
----------
nPts: int
Number of points to be sampled
spacingFunc: function
sampling function object. The methods available in :mod:`prefoil.sampling` are a good default example.
func_args: dictionary
Dictionary of input arguments for the sampling function.
nTEPts: float
Number of points along the **blunt** trailing edge
TE_knot: bool
If True, add a duplicate point between the lower airfoil surface and the TE to indicate that a knot is present. If there is a sharp or round trailing edge then this does nothing.
Returns
-------
coords : Ndarray [N, 2]
Coordinates array, anticlockwise, from trailing edge
"""
s = sampling.joinedSpacing(nPts, spacingFunc=spacingFunc, func_args=func_args, s_LE=self.s_LE)
sampled_coords = self.spline.getValue(s)
if not self.closedCurve and TE_knot:
sampled_coords = np.vstack((sampled_coords, sampled_coords[-1]))
if nTEPts and not self.closedCurve:
coords_TE = np.zeros((nTEPts + 2, sampled_coords.shape[1]))
for idim in range(sampled_coords.shape[1]):
val1 = self.spline.getValue(1)[idim]
val2 = self.spline.getValue(0)[idim]
coords_TE[:, idim] = np.linspace(val1, val2, nTEPts + 2)
sampled_coords = np.vstack((sampled_coords, coords_TE[1:-1]))
if not self.closedCurve:
sampled_coords = np.vstack((sampled_coords, sampled_coords[0]))
self.sampled_pts = sampled_coords
return sampled_coords
def _buildFFD(self, nffd, fitted, xmargin, ymarginu, ymarginl, xslice, coords):
"""
The function that actually builds the FFD Box from all of the given parameters
Parameters
----------
nffd : int
number of FFD points along the chord
fitted : bool
flag to pick between a fitted FFD (True) and box FFD (False)
xmargin : float
The closest distance of the FFD box to the tip and aft of the airfoil
ymarginu : float
When a box ffd is generated this specifies the top of the box's y values as
the maximum y value in the airfoil coordinates plus this margin.
When a fitted ffd is generated this is the margin between the FFD point at
an xslice location and the upper surface of the airfoil at this location
ymarginl : float
When a box ffd is generated this specifies the bottom of the box's y values as
the minimum y value in the airfoil coordinates minus this margin.
When a fitted ffd is generated this is the margin between the FFD point at
an xslice location and the lower surface of the airfoil at this location
xslice : Ndarray [N,2]
User specified xslice locations. If this is chosen nffd is ignored
coords : Ndarray [N,2]
the coordinates to use for defining the airfoil, if the user does not
want the original coordinates for the airfoil used. This shouldn't be
used unless the user wants fine tuned control over the FFD creation,
It should be sufficient to ignore.
"""
if coords is None:
coords = self.getSplinePts()
if xslice is None:
xslice = np.zeros(nffd)
for i in range(nffd):
xtemp = i * 1.0 / (nffd - 1.0)
xslice[i] = min(coords[:, 0]) - 1.0 * xmargin + (max(coords[:, 0]) + 2.0 * xmargin) * xtemp
else:
nffd = len(xslice)
FFDbox = np.zeros((nffd, 2, 2, 3))
if fitted:
ylower = np.zeros(nffd)
yupper = np.zeros(nffd)
for i in range(nffd):
ymargin = ymarginu + (ymarginl - ymarginu) * xslice[i]
yu, yl = _getClosestY(coords, xslice[i])
yupper[i] = yu + ymargin
ylower[i] = yl - ymargin
else:
yupper = np.ones(nffd) * (max(coords[:, 1]) + ymarginu)
ylower = np.ones(nffd) * (min(coords[:, 1]) - ymarginl)
# X
FFDbox[:, 0, 0, 0] = xslice[:].copy()
FFDbox[:, 1, 0, 0] = xslice[:].copy()
# Y
# lower
FFDbox[:, 0, 0, 1] = ylower[:].copy()
# upper
FFDbox[:, 1, 0, 1] = yupper[:].copy()
# copy
FFDbox[:, :, 1, :] = FFDbox[:, :, 0, :].copy()
# Z
FFDbox[:, :, 0, 2] = 0.0
# Z
FFDbox[:, :, 1, 2] = 1.0
return FFDbox
## Output
[docs] def writeCoords(self, filename, coords=None, spline_coords=False, file_format="plot3d"):
"""
Writes out a set of airfoil coordinates. By default, the most recently sampled coordinates are written out.
If there are no recently sampled coordinates and none are passed in this will fail.
Parameters
----------
filename : str
the filename without extension to write to
coords : Ndarray [N,2]
the coordinates to write out to a file. If None then the most recent sampled set of points is used.
spline_coords : bool
If true it will write out the underlying spline coordinates and the value of `coords` will be ignored. Useful if only geometric modifications to coordinates are being preformed.
file_format : str
the file format to write, can be `plot3d` or `dat`
"""
if spline_coords is True:
coords = self.getSplinePts()
if coords is None:
if self.sampled_pts is not None:
coords = self.sampled_pts
else:
raise Error("No coordinates to write!")
if file_format == "plot3d":
_writePlot3D(filename, coords[:, 0], coords[:, 1])
elif file_format == "dat":
_writeDat(filename, coords[:, 0], coords[:, 1])
else:
raise Error(file_format + " is not a supported output format!")
[docs] def generateFFD(
self, nffd, filename, fitted=True, xmargin=0.001, ymarginu=0.02, ymarginl=0.02, xslice=None, coords=None
):
"""
Generates an FFD from the airfoil and writes it out to file
Parameters
----------
nffd : int
the number of chordwise points in the FFD
filename : str
filename to write out, not including the '.xyz' ending
fitted : bool
flag to pick between a fitted FFD (True) and box FFD (False)
xmargin : float
The closest distance of the FFD box to the tip and aft of the airfoil
ymarginu : float
When a box ffd is generated this specifies the top of the box's y values as
the maximum y value in the airfoil coordinates plus this margin.
When a fitted ffd is generated this is the margin between the FFD point at
an xslice location and the upper surface of the airfoil at this location
ymarginl : float
When a box ffd is generated this specifies the bottom of the box's y values as
the minimum y value in the airfoil coordinates minus this margin.
When a fitted ffd is generated this is the margin between the FFD point at
an xslice location and the lower surface of the airfoil at this location
xslice : Ndarray [N,2]
User specified xslice locations. If this is chosen nffd is ignored
coords : Ndarray [N,2]
the coordinates to use for defining the airfoil, if the user does not
want the original coordinates for the airfoil used. This shouldn't be
used unless the user wants fine tuned control over the FFD creation,
It should be sufficient to ignore.
"""
FFDbox = self._buildFFD(nffd, fitted, xmargin, ymarginu, ymarginl, xslice, coords)
_writeFFD(FFDbox, filename)
## Utils
# maybe remove and put into a separate location?
[docs] def plot(self, camber=False):
"""
Plots the airfoil.
It tries to plot the most recently sampled set of points, but if none exists, it will plot the original set of coordinates.
Parameters
----------
camber : bool
True to plot the camber line
Returns
-------
fig : matplotlib.pyplot.Figure
The figure with the plotted airfoil
"""
try:
import matplotlib.pyplot as plt
except ImportError as e:
raise ImportError("matplotlib is needed for the plotting functionality") from e
if self.sampled_pts is None:
coords = self.getSplinePts()
else:
coords = self.sampled_pts
fig = plt.figure()
# pts = self._getDefaultSampling(npts=1000)
plt.plot(coords[:, 0], coords[:, 1], "-r")
plt.axis("equal")
# if self.sampled_X is not None:
plt.plot(coords[:, 0], coords[:, 1], "o")
if camber:
camber_pts = self.camber.getValue(np.linspace(0, 1, 200))
plt.plot(camber_pts[:, 0], camber_pts[:, 1], "--g", label="camber")
return fig