diff --git a/bladex/profilebase.py b/bladex/profilebase.py index a0a1917..fb0eecc 100644 --- a/bladex/profilebase.py +++ b/bladex/profilebase.py @@ -1,3 +1,7 @@ +""" +Base module for blade construction. Provide essential tools for on airfoils. +""" + import numpy as np import matplotlib.pyplot as plt from .ndinterpolator import reconstruct_f @@ -7,15 +11,30 @@ 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 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 + 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): @@ -23,17 +42,29 @@ def __init__(self): 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) 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. + 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 self.xup_coordinates[0] != self.xdown_coordinates[0]: + raise ValueError('Airfoils must have xup_coordinates[0] \ + != xdown_coordinates[0]') + if self.xup_coordinates[-1] != self.xdown_coordinates[-1]: + raise ValueError('Airfoils must have xup_coordinates[-1] \ + != 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] @@ -44,52 +75,34 @@ def _update_edges(self): self.trailing_edge[1] = 0.5 * ( self.yup_coordinates[-1] + self.ydown_coordinates[-1]) - @property - def reference_point(self): - """ - Returns the coordinates of the airfoil's geometric center. - - :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): - """ - Measures 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) - def interpolate_coordinates(self, num=500, radius=1.0): """ - Interpolates the airfoil coordinates from the given data set of discrete points. + 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 :doc:`/utils/rbf`. + 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 :doc:`/utils/rbf`. References: - Buhmann, Martin D. (2003), Radial Basis Functions: Theory and Implementations. + 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 + :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 + :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.') @@ -119,82 +132,280 @@ def interpolate_coordinates(self, num=500, radius=1.0): return xx_up, xx_down, yy_up, yy_down - def max_thickness(self, interpolate=False): + def compute_chord_line(self, interpolate=False, n_interpolated_points=500): + """ + 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 bool interpolate: if True, then the coordinates are computed at + equally spaced intervals within the range of LE and TE. Default + value is False + :param int n_interpolated_points: number of points to be used for the + equally-spaced sample computations, and is used only if parameter + interpolate is True. Default value is 500 """ - Returns the airfoil's maximum thickness. + self._update_edges() + aratio = ((self.trailing_edge[1] - self.leading_edge[1]) / + (self.trailing_edge[0] - self.leading_edge[0])) + if ((interpolate is True) or + (self.xup_coordinates == self.xdown_coordinates).all() is False): + cl_x_coordinates = np.linspace( + self.leading_edge[0], + self.trailing_edge[0], + num=n_interpolated_points) + cl_y_coordinates = np.zeros(cl_x_coordinates.size) + 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 = np.zeros(self.xup_coordinates.size) + 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]) - 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. + def compute_camber_line(self, interpolate=False, n_interpolated_points=500): + """ + 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 bool interpolate: if True, the interpolated coordinates are used + to obtain the camber line, otherwise the original discrete + coordinates are used. Default value is False. + :param int n_interpolated_points: number of points to be used for the + uniform interpolation, and is used only if parameter interpolate + is True. Default value is 500 + + 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 camberline. + """ + if (interpolate is True) or ( + (self.xup_coordinates == self.xdown_coordinates).all() is False): + # 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. + 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]) + + def deform_camber_line(self, + percent_change, + interpolate=False, + n_interpolated_points=500): + """ + Deform camber line according to a given percentage of change of the + maximum camber. Also reconstructs the deformed airfoil's coordinates. - 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. + The percentage of change is defined as follows: + .. math:: + \\frac{new magnitude of max camber - old magnitude of maximum \ + camber}{old magnitude of maximum camber} * 100 + + A positive percentage means the new camber is larger than the max + camber value, while a negative percentage indicates the new value + is smaller. + + We note that the method works only for airfoils in the reference + position, i.e. chord line lies on the X-axis and the foil is not + rotated, since the measurements are based on the Y-values of the + airfoil coordinates, hence any measurements or scalings will be + inaccurate for the foils not in their reference position. + + :param float max_camber_change_percent: percentage of change of the + maximum camber. Default value is None + :param bool interpolate: if True, the interpolated coordinates are + used to compute the camber line and foil's thickness, otherwise + the original discrete coordinates are used. Default value is False. + :param n_interpolated_points: number of points to be used for the + uniform interpolation, and is used only if parameter interpolate + is True. Default value is 500 + :raises ValueError: if max_camber_change_percent is None + """ + # Updating camber line + self.compute_camber_line( + interpolate=interpolate, + n_interpolated_points=n_interpolated_points) + scaling_factor = percent_change / 100. + 1. + self.camber_line[1] *= scaling_factor + + # Evaluating half-thickness of the undeformed airfoil, + # which should hold same values for the deformed foil. + if (interpolate is True) or ( + (self.xup_coordinates == self.xdown_coordinates).all() is False): + # 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. + (self.xup_coordinates, self.xdown_coordinates, self.yup_coordinates, + self.ydown_coordinates + ) = self.interpolate_coordinates(num=n_interpolated_points) + + half_thickness = 0.5 * np.fabs( + self.yup_coordinates - self.ydown_coordinates) + + self.yup_coordinates = self.camber_line[1] + half_thickness + self.ydown_coordinates = self.camber_line[1] - half_thickness - :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 + @property + def reference_point(self): + """ + Return the coordinates of the airfoil's geometric center. + + :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) + + def get_max_thickness(self, interpolate=False, n_interpolated_points=500): + """ + 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 + uniform interpolation, and is used only if parameter interpolate is + True. Default value is 500 :return: maximum thickness :rtype: float """ if (interpolate is True) or ( (self.xup_coordinates == self.xdown_coordinates).all() is False): - # Evaluation of the thickness requires comparing both y_up and y_down for the same x-section, - # (i.e. same x_coordinate), according to the 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. - xx_up, xx_down, yy_up, yy_down = self.interpolate_coordinates() + # 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() - def max_camber(self, interpolate=False): + def get_max_camber(self, interpolate=False, n_interpolated_points=500): """ - Returns 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 + 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 + uniform interpolation, and is used only if parameter interpolate + is True. Default value is 500 :return: maximum camber :rtype: float """ - if (interpolate is True) or ( - (self.xup_coordinates == self.xdown_coordinates).all() is False): - # Evaluation of camber requires comparing both y_up and y_down for the same x-section (i.e. same x_coordinate). - # 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. - xx_up, xx_down, yy_up, yy_down = self.interpolate_coordinates() - camber = yy_down + 0.5 * np.fabs(yy_up - yy_down) + self.compute_chord_line( + interpolate=interpolate, + n_interpolated_points=n_interpolated_points) + + self.compute_camber_line( + interpolate=interpolate, + 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]) + + if (self.camber_line[1][camber.argmax()] > + self.chord_line[1][camber.argmax()]): + # Camber line is above the chord line, at the point of max camber + max_camber = camber.max() else: - camber = self.ydown_coordinates + 0.5 * np.fabs( - self.yup_coordinates - self.ydown_coordinates) - return camber.max() + # Camber line is below the chord line, at the point of max camber + max_camber = camber.max() * -1 + + return max_camber def rotate(self, rad_angle=None, deg_angle=None): """ - 2D counter clockwise rotation about the origin of the Cartesian coordinate system. + 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. + 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) = \\left(\\begin{matrix} x \\ y \\end{matrix}\\right)`, the rotated coordinates will be: - .. math:: - P^{'} = \\left(\\begin{matrix} x^{'} \\ y^{'} \\end{matrix}\\right) = R (\\theta) \\cdot P` + Given the coordinates of point :math:`(P) = + \\left(\\begin{matrix} x \\ y \\end{matrix}\\right)`, + 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. + 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 + :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.') + 'You have to pass either the angle in radians or in degrees, \ + not both.') if rad_angle: cosine = np.cos(rad_angle) sine = np.sin(rad_angle) @@ -228,7 +439,7 @@ def rotate(self, rad_angle=None, deg_angle=None): def translate(self, translation): """ - Translates the airfoil coordinates according to a 2D translation vector. + Translate the airfoil coordinates according to a 2D translation vector. :param array_like translation: the translation vector in 2D """ @@ -239,7 +450,7 @@ def translate(self, translation): def flip(self): """ - Flips the airfoil coordinates about both the X-axis and the Y-axis. + Flip the airfoil coordinates about both the X-axis and the Y-axis. """ self.xup_coordinates *= -1 self.xdown_coordinates *= -1 @@ -248,10 +459,13 @@ def flip(self): def scale(self, factor): """ - Scales the airfoil coordinates according to a scaling 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. + 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 """ @@ -265,9 +479,11 @@ def scale(self, factor): def plot(self, outfile=None): """ - Plots the airfoil coordinates. + Plot the airfoil coordinates. - :param string 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 + :param string 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() plt.plot(self.xup_coordinates, self.yup_coordinates) @@ -276,6 +492,6 @@ def plot(self, outfile=None): plt.axis('equal') if outfile: - if not isinstance(outfile, string): + if not isinstance(outfile, str): raise ValueError('Output file name must be string.') plt.savefig(outfile)