"""
Base module that provides essential tools and transformations on airfoils.
"""
import numpy as np
import matplotlib.pyplot as plt
from .ndinterpolator import reconstruct_f
[docs]class ProfileBase(object):
"""
Base sectional profile of the propeller blade.
Each sectional profile is a 2D airfoil that is split into two parts: the
upper and lower parts. The coordiates of each part is represented by two
arrays corresponding to the X and Y components in the 2D coordinate system.
Such coordinates can be either generated using NACA functions, or be
inserted directly by the user as custom profiles.
:param numpy.ndarray xup_coordinates: 1D array that contains the
X-components of the airfoil upper-half surface. Default value is None
:param numpy.ndarray xdown_coordinates: 1D array that contains the
X-components of the airfoil lower-half surface. Default value is None
:param numpy.ndarray yup_coordinates: 1D array that contains the
Y-components of the airfoil upper-half surface. Default value is None
:param numpy.ndarray ydown_coordinates: 1D array that contains the
Y-components of the airfoil lower-half surface. Default value is None
:param numpy.ndarray chord_line: contains the X and Y coordinates of the
straight line joining between the leading and trailing edges. Default
value is None
:param numpy.ndarray camber_line: contains the X and Y coordinates of the
curve passing through all the mid-points between the upper and lower
surfaces of the airfoil. Default value is None
:param numpy.ndarray leading_edge: 2D coordinates of the airfoil's
leading edge. Default values are zeros
:param numpy.ndarray trailing_edge: 2D coordinates of the airfoil's
trailing edge. Default values are zeros
"""
def __init__(self):
self.xup_coordinates = None
self.xdown_coordinates = None
self.yup_coordinates = None
self.ydown_coordinates = None
self.chord_line = None
self.camber_line = None
self.leading_edge = np.zeros(2)
self.trailing_edge = np.zeros(2)
[docs] def _update_edges(self):
"""
Private method that identifies and updates the airfoil's leading and
trailing edges.
Given the airfoil coordinates from the leading to the trailing edge,
if the trailing edge has a non-zero thickness, then the average value
between the upper and lower trailing edges is taken as the true
trailing edge, hence both the leading and the trailing edges are always
unique.
"""
if np.fabs(self.xup_coordinates[0] - self.xdown_coordinates[0]) > 1e-4:
raise ValueError('Airfoils must have xup_coordinates[0] '\
'almost equal to xdown_coordinates[0]')
if np.fabs(
self.xup_coordinates[-1] - self.xdown_coordinates[-1]) > 1e-4:
raise ValueError('Airfoils must have xup_coordinates[-1] '\
'almost equal to xdown_coordinates[-1]')
self.leading_edge[0] = self.xup_coordinates[0]
self.leading_edge[1] = self.yup_coordinates[0]
self.trailing_edge[0] = self.xup_coordinates[-1]
if self.yup_coordinates[-1] == self.ydown_coordinates[-1]:
self.trailing_edge[1] = self.yup_coordinates[-1]
else:
self.trailing_edge[1] = 0.5 * (
self.yup_coordinates[-1] + self.ydown_coordinates[-1])
[docs] def interpolate_coordinates(self, num=500, radius=1.0):
"""
Interpolate the airfoil coordinates from the given data set of
discrete points.
The interpolation applies the Radial Basis Function (RBF) method,
to construct approximations of the two functions that correspond to the
airfoil upper half and lower half coordinates. The RBF implementation
is present in :ref:`RBF ndinterpolator <ndinterpolator-label>`.
References:
Buhmann, Martin D. (2003), Radial Basis Functions: Theory
and Implementations.
http://www.cs.bham.ac.uk/~jxb/NN/l12.pdf
https://www.cc.gatech.edu/~isbell/tutorials/rbf-intro.pdf
:param int num: number of interpolated points. Default value is 500
:param float radius: range of the cut-off radius necessary for the RBF
interpolation. Default value is 1.0. It is quite necessary to
adjust the value properly so as to ensure a smooth interpolation
:return: interpolation points for the airfoil upper half X-component,
interpolation points for the airfoil lower half X-component,
interpolation points for the airfoil upper half Y-component,
interpolation points for the airfoil lower half Y-component
:rtype: numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray
:raises TypeError: if num is not of type int
:raises ValueError: if num is not positive, or if radius is not
positive
"""
if not isinstance(num, int):
raise TypeError('Inserted value must be of type integer.')
if num <= 0 or radius <= 0:
raise ValueError('Inserted value must be positive.')
xx_up = np.linspace(
self.xup_coordinates[0], self.xup_coordinates[-1], num=num)
yy_up = np.zeros(num)
reconstruct_f(
basis='beckert_wendland_c2_basis',
radius=radius,
original_input=self.xup_coordinates,
original_output=self.yup_coordinates,
rbf_input=xx_up,
rbf_output=yy_up)
xx_down = np.linspace(
self.xdown_coordinates[0], self.xdown_coordinates[-1], num=num)
yy_down = np.zeros(num)
reconstruct_f(
basis='beckert_wendland_c2_basis',
radius=radius,
original_input=self.xdown_coordinates,
original_output=self.ydown_coordinates,
rbf_input=xx_down,
rbf_output=yy_down)
return xx_up, xx_down, yy_up, yy_down
[docs] def compute_chord_line(self, n_interpolated_points=None):
"""
Compute the 2D coordinates of the chord line. Also updates
the chord_line class member.
The chord line is the straight line that joins between the leading edge
and the trailing edge. It is simply computed from the equation of
a line passing through two points, the LE and TE.
:param int n_interpolated_points: number of points to be used for the
equally-spaced sample computations. If None then there is no
interpolation, unless the arrays x_up != x_down elementwise which
implies that the corresponding y_up and y_down can not be
comparable, hence a uniform interpolation is required. Default
value is None
"""
self._update_edges()
aratio = ((self.trailing_edge[1] - self.leading_edge[1]) /
(self.trailing_edge[0] - self.leading_edge[0]))
if not (self.xup_coordinates == self.xdown_coordinates
).all() and n_interpolated_points is None:
# If x_up != x_down element-wise, then the corresponding y_up and
# y_down can not be comparable, hence a uniform interpolation is
# required. Also in case the interpolated_points is None,
# then we assume a default number of interpolated points
n_interpolated_points = 500
if n_interpolated_points:
cl_x_coordinates = np.linspace(
self.leading_edge[0],
self.trailing_edge[0],
num=n_interpolated_points)
cl_y_coordinates = (aratio *
(cl_x_coordinates - self.leading_edge[0]) +
self.leading_edge[1])
self.chord_line = np.array([cl_x_coordinates, cl_y_coordinates])
else:
cl_y_coordinates = (aratio *
(self.xup_coordinates - self.leading_edge[0]) +
self.leading_edge[1])
self.chord_line = np.array([self.xup_coordinates, cl_y_coordinates])
[docs] def compute_camber_line(self, n_interpolated_points=None):
"""
Compute the 2D coordinates of the camber line. Also updates the
camber_line class member.
The camber line is defined by the curve passing through all the mid
points between the upper surface and the lower surface of the airfoil.
:param int n_interpolated_points: number of points to be used for the
equally-spaced sample computations. If None then there is no
interpolation, unless the arrays x_up != x_down elementwise which
implies that the corresponding y_up and y_down can not be
comparable, hence a uniform interpolation is required. Default
value is None
We note that a uniform interpolation becomes necessary for the cases
when the X-coordinates of the upper and lower surfaces do not
correspond to the same vertical sections, since this would imply
inaccurate measurements for obtaining the camber line.
"""
if not (self.xup_coordinates == self.xdown_coordinates
).all() and n_interpolated_points is None:
# If x_up != x_down element-wise, then the corresponding y_up and
# y_down can not be comparable, hence a uniform interpolation is
# required. Also in case the interpolated_points is None,
# then we assume a default number of interpolated points
n_interpolated_points = 500
if n_interpolated_points:
cl_x_coordinates, yy_up, yy_down = (
self.interpolate_coordinates(num=n_interpolated_points)[1:])
cl_y_coordinates = 0.5 * (yy_up + yy_down)
self.camber_line = np.array([cl_x_coordinates, cl_y_coordinates])
else:
cl_y_coordinates = (0.5 *
(self.ydown_coordinates + self.yup_coordinates))
self.camber_line = np.array(
[self.xup_coordinates, cl_y_coordinates])
@property
def reference_point(self):
"""
Return the coordinates of the chord's mid point.
:return: reference point in 2D
:rtype: numpy.ndarray
"""
self._update_edges()
reference_point = [
0.5 * (self.leading_edge[0] + self.trailing_edge[0]),
0.5 * (self.leading_edge[1] + self.trailing_edge[1])
]
return np.asarray(reference_point)
@property
def chord_length(self):
"""
Measure the l2-norm (Euclidean distance) between the leading edge
and the trailing edge.
:return: chord length
:rtype: float
"""
self._update_edges()
return np.linalg.norm(self.leading_edge - self.trailing_edge)
[docs] def max_thickness(self, n_interpolated_points=None):
"""
Return the airfoil's maximum thickness.
Thickness is defined as the distnace between the upper and lower
surfaces of the airfoil, and can be measured in two different ways:
- American convention: measures along the line perpendicular to \
the mean camber line.
- British convention: measures along the line perpendicular to \
the chord line.
In this implementation, the british convention is used to evaluate
the maximum thickness.
References:
Phillips, Warren F. (2010). Mechanics of Flight (2nd ed.). \
Wiley & Sons. p. 27. ISBN 978-0-470-53975-0.
Bertin, John J.; Cummings, Russel M. (2009). Pearson Prentice Hall, \
ed. Aerodynamics for Engineers (5th ed.). \
p. 199. ISBN 978-0-13-227268-1.
:param bool interpolate: if True, the interpolated coordinates are used
to measure the thickness; otherwise, the original discrete
coordinates are used. Default value is False
:param int n_interpolated_points: number of points to be used for the
equally-spaced sample computations. If None then there is no
interpolation, unless the arrays x_up != x_down elementwise which
implies that the corresponding y_up and y_down can not be
comparable, hence a uniform interpolation is required. Default
value is None
:return: maximum thickness
:rtype: float
"""
if not (self.xup_coordinates == self.xdown_coordinates
).all() and n_interpolated_points is None:
# If x_up != x_down element-wise, then the corresponding y_up and
# y_down can not be comparable, hence a uniform interpolation is
# required. Also in case the interpolated_points is None,
# then we assume a default number of interpolated points
n_interpolated_points = 500
if n_interpolated_points:
# Evaluation of the thickness requires comparing both y_up and
# y_down for the same x-section, (i.e. same x_coordinate),
# according to british convention. If x_up != x_down element-wise,
# then the corresponding y_up and y_down can not be comparable,
# hence a uniform interpolation is required.
yy_up, yy_down = self.interpolate_coordinates(
num=n_interpolated_points)[2:]
return np.fabs(yy_up - yy_down).max()
return np.fabs(self.yup_coordinates - self.ydown_coordinates).max()
[docs] def max_camber(self, n_interpolated_points=500):
"""
Return the magnitude of the airfoil's maximum camber.
Camber is defined as the distance between the chord line and the mean
camber line, and is measured along the line perpendicular to the chord
line.
:param bool interpolate: if True, the interpolated coordinates are used
to measure the camber; otherwise, the original discrete coordinates
are used. Default value is False
:param int n_interpolated_points: number of points to be used for the
equally-spaced sample computations. If None then there is no
interpolation, unless the arrays x_up != x_down elementwise which
implies that the corresponding y_up and y_down can not be
comparable, hence a uniform interpolation is required. Default
value is None
:return: maximum camber
:rtype: float
"""
self.compute_chord_line(n_interpolated_points=n_interpolated_points)
self.compute_camber_line(n_interpolated_points=n_interpolated_points)
n_points = self.camber_line[0].size
camber = np.zeros(n_points)
for i in range(n_points):
camber[i] = np.linalg.norm(
self.chord_line[:, i] - self.camber_line[:, i])
max_camber = camber.max()
if (self.camber_line[1][camber.argmax()] <
self.chord_line[1][camber.argmax()]):
# Camber line is below the chord line, at the point of max camber
max_camber *= -1
return max_camber
[docs] def rotate(self, rad_angle=None, deg_angle=None):
"""
2D counter clockwise rotation about the origin of the Cartesian
coordinate system.
The rotation matrix, :math:`R(\\theta)`, is used to perform rotation
in the 2D Euclidean space about the origin, which is -- by default --
the leading edge.
:math:`R(\\theta)` is defined by:
.. math::
\\left(\\begin{matrix} cos (\\theta) & - sin (\\theta) \\\\
sin (\\theta) & cos (\\theta) \\end{matrix}\\right)
Given the coordinates of point :math:`P` such that
.. math::
P = \\left(\\begin{matrix} x \\\\
y \\end{matrix}\\right),
Then, the rotated coordinates will be:
.. math::
P^{'} = \\left(\\begin{matrix} x^{'} \\\\
y^{'} \\end{matrix}\\right)
= R (\\theta) \\cdot P
If a standard right-handed Cartesian coordinate system is used, with
the X-axis to the right and the Y-axis up, the rotation
:math:`R (\\theta)` is counterclockwise. If a left-handed Cartesian
coordinate system is used, with X-axis directed to the right and Y-axis
directed down, :math:`R (\\theta)` is clockwise.
:param float rad_angle: angle in radians. Default value is None
:param float deg_angle: angle in degrees. Default value is None
:raises ValueError: if both rad_angle and deg_angle are inserted,
or if neither is inserted
"""
if rad_angle is not None and deg_angle is not None:
raise ValueError(
'You have to pass either the angle in radians or in degrees,' \
' not both.')
if rad_angle is not None:
cosine = np.cos(rad_angle)
sine = np.sin(rad_angle)
elif deg_angle is not None:
cosine = np.cos(np.radians(deg_angle))
sine = np.sin(np.radians(deg_angle))
else:
raise ValueError(
'You have to pass either the angle in radians or in degrees.')
rot_matrix = np.array([cosine, -sine, sine, cosine]).reshape((2, 2))
coord_matrix_up = np.vstack((self.xup_coordinates,
self.yup_coordinates))
coord_matrix_down = np.vstack((self.xdown_coordinates,
self.ydown_coordinates))
new_coord_matrix_up = np.zeros(coord_matrix_up.shape)
new_coord_matrix_down = np.zeros(coord_matrix_down.shape)
for i in range(self.xup_coordinates.shape[0]):
new_coord_matrix_up[:, i] = np.dot(rot_matrix,
coord_matrix_up[:, i])
new_coord_matrix_down[:, i] = np.dot(rot_matrix,
coord_matrix_down[:, i])
self.xup_coordinates = new_coord_matrix_up[0]
self.xdown_coordinates = new_coord_matrix_down[0]
self.yup_coordinates = new_coord_matrix_up[1]
self.ydown_coordinates = new_coord_matrix_down[1]
[docs] def translate(self, translation):
"""
Translate the airfoil coordinates according to a 2D translation vector.
:param array_like translation: the translation vector in 2D
"""
self.xup_coordinates += translation[0]
self.xdown_coordinates += translation[0]
self.yup_coordinates += translation[1]
self.ydown_coordinates += translation[1]
[docs] def reflect(self):
"""
Reflect the airfoil coordinates about the origin, i.e. a mirror
transformation is performed about both the X-axis and the Y-axis.
"""
self.xup_coordinates *= -1
self.xdown_coordinates *= -1
self.yup_coordinates *= -1
self.ydown_coordinates *= -1
[docs] def scale(self, factor):
"""
Scale the airfoil coordinates according to a scaling factor.
In order to apply the scaling without affecting the position of the
reference point, the method translates the airfoil by its refernce
point to be centered in the origin, then the scaling is applied, and
finally the airfoil is translated back by its reference point to the
initial position.
:param float factor: the scaling factor
"""
ref_point = self.reference_point
self.translate(-ref_point)
self.xup_coordinates *= factor
self.xdown_coordinates *= factor
self.yup_coordinates *= factor
self.ydown_coordinates *= factor
self.translate(ref_point)
[docs] def plot(self,
profile=True,
chord_line=False,
camber_line=False,
ref_point=False,
outfile=None):
"""
Plot the airfoil coordinates.
:param bool profile: if True, then plot the profile coordinates.
Default value is True
:param bool chord_line: if True, then plot the chord line. Default
value is False
:param bool camber_line: if True, then plot the camber line. Default
value is False
:param bool ref_point: if True, then scatter plot the reference point.
Default value is False
:param str outfile: outfile name. If a string is provided then the
plot is saved with that name, otherwise the plot is not saved.
Default value is None
"""
plt.figure()
if (self.xup_coordinates is None or self.yup_coordinates is None
or self.xdown_coordinates is None
or self.ydown_coordinates is None):
raise ValueError('One or all the coordinates have None value.')
if profile:
plt.plot(
self.xup_coordinates,
self.yup_coordinates,
label='Upper profile')
plt.plot(
self.xdown_coordinates,
self.ydown_coordinates,
label='Lower profile')
if chord_line:
if self.chord_line is None:
raise ValueError(
'Chord line is None. You must compute it first')
plt.plot(self.chord_line[0], self.chord_line[1], label='Chord line')
if camber_line:
if self.camber_line is None:
raise ValueError(
'Camber line is None. You must compute it first')
plt.plot(
self.camber_line[0], self.camber_line[1], label='Camber line')
if ref_point:
plt.scatter(
self.reference_point[0],
self.reference_point[1],
s=15,
label='Reference point')
plt.grid(linestyle='dotted')
plt.axis('equal')
plt.legend()
if outfile:
if not isinstance(outfile, str):
raise ValueError('Output file name must be string.')
plt.savefig(outfile)
else:
plt.show()