diff --git a/.pylintrc b/.pylintrc index 80140b391..dbb2cb11d 100644 --- a/.pylintrc +++ b/.pylintrc @@ -215,7 +215,7 @@ good-names=FlightPhases, prop_I_11, Kt, # transformation matrix transposed clalpha2D, - clalpha2D_incompresible, + clalpha2D_incompressible, r_NOZ, # Nozzle position vector rocket_dry_I_33, rocket_dry_I_11, diff --git a/CHANGELOG.md b/CHANGELOG.md index 6cbf037aa..eb2bfec8f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,7 @@ Attention: The newest changes should be on top --> ### Changed - MNT: Refactors the code to adopt pylint [#621](https://github.com/RocketPy-Team/RocketPy/pull/621) +- MNT: Refactor AeroSurfaces [#634](https://github.com/RocketPy-Team/RocketPy/pull/634) ### Fixed diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py deleted file mode 100644 index 582fdb660..000000000 --- a/rocketpy/rocket/aero_surface.py +++ /dev/null @@ -1,2149 +0,0 @@ -import warnings -from abc import ABC, abstractmethod - -import numpy as np -from scipy.optimize import fsolve - -from ..mathutils.function import Function -from ..plots.aero_surface_plots import ( - _AirBrakesPlots, - _EllipticalFinsPlots, - _NoseConePlots, - _TailPlots, - _TrapezoidalFinsPlots, -) -from ..prints.aero_surface_prints import ( - _AirBrakesPrints, - _EllipticalFinsPrints, - _NoseConePrints, - _RailButtonsPrints, - _TailPrints, - _TrapezoidalFinsPrints, -) - -# TODO: all the evaluate_shape() methods need tests and documentation - - -class AeroSurface(ABC): - """Abstract class used to define aerodynamic surfaces.""" - - def __init__(self, name): - self.cpx = 0 - self.cpy = 0 - self.cpz = 0 - self.name = name - - @staticmethod - def _beta(mach): - """Defines a parameter that is often used in aerodynamic - equations. It is commonly used in the Prandtl factor which - corrects subsonic force coefficients for compressible flow. - This is applied to the lift coefficient of the nose cone, - fins and tails/transitions as in [1]. - - Parameters - ---------- - mach : int, float - Number of mach. - - Returns - ------- - beta : int, float - Value that characterizes flow speed based on the mach number. - - References - ---------- - [1] Barrowman, James S. https://arc.aiaa.org/doi/10.2514/6.1979-504 - """ - - if mach < 0.8: - return np.sqrt(1 - mach**2) - elif mach < 1.1: - return np.sqrt(1 - 0.8**2) - else: - return np.sqrt(mach**2 - 1) - - @abstractmethod - def evaluate_center_of_pressure(self): - """Evaluates the center of pressure of the aerodynamic surface in local - coordinates. - - Returns - ------- - None - """ - - @abstractmethod - def evaluate_lift_coefficient(self): - """Evaluates the lift coefficient curve of the aerodynamic surface. - - Returns - ------- - None - """ - - @abstractmethod - def evaluate_geometrical_parameters(self): - """Evaluates the geometrical parameters of the aerodynamic surface. - - Returns - ------- - None - """ - - @abstractmethod - def info(self): - """Prints and plots summarized information of the aerodynamic surface. - - Returns - ------- - None - """ - - @abstractmethod - def all_info(self): - """Prints and plots all the available information of the aero surface. - - Returns - ------- - None - """ - - -class NoseCone(AeroSurface): - """Keeps nose cone information. - - Note - ---- - The local coordinate system has the origin at the tip of the nose cone - and the Z axis along the longitudinal axis of symmetry, positive - downwards (top -> bottom). - - Attributes - ---------- - NoseCone.length : float - Nose cone length. Has units of length and must be given in meters. - NoseCone.kind : string - Nose cone kind. Can be "conical", "ogive", "elliptical", "tangent", - "von karman", "parabolic", "powerseries" or "lvhaack". - NoseCone.bluffness : float - Ratio between the radius of the circle on the tip of the ogive and the - radius of the base of the ogive. Currently only used for the nose cone's - drawing. Must be between 0 and 1. Default is None, which means that the - nose cone will not have a sphere on the tip. If a value is given, the - nose cone's length will be slightly reduced because of the addition of - the sphere. Must be None or 0 if a "powerseries" nose cone kind is - specified. - NoseCone.rocket_radius : float - The reference rocket radius used for lift coefficient normalization, - in meters. - NoseCone.base_radius : float - Nose cone base radius. Has units of length and must be given in meters. - NoseCone.radius_ratio : float - Ratio between the nose cone base radius and the rocket radius. Has no - units. If base radius is not given, the ratio between base radius and - rocket radius is assumed as 1, meaning that the nose cone has the same - radius as the rocket. If base radius is given, the ratio between base - radius and rocket radius is calculated and used for lift calculation. - NoseCone.power : float - Factor that controls the bluntness of the shape for a power series - nose cone. Must be between 0 and 1. It is ignored when other nose - cone types are used. - NoseCone.name : string - Nose cone name. Has no impact in simulation, as it is only used to - display data in a more organized matter. - NoseCone.cp : tuple - Tuple with the x, y and z local coordinates of the nose cone center of - pressure. Has units of length and is given in meters. - NoseCone.cpx : float - Nose cone local center of pressure x coordinate. Has units of length and - is given in meters. - NoseCone.cpy : float - Nose cone local center of pressure y coordinate. Has units of length and - is given in meters. - NoseCone.cpz : float - Nose cone local center of pressure z coordinate. Has units of length and - is given in meters. - NoseCone.cl : Function - Function which defines the lift coefficient as a function of the angle - of attack and the Mach number. Takes as input the angle of attack in - radians and the Mach number. Returns the lift coefficient. - NoseCone.clalpha : float - Lift coefficient slope. Has units of 1/rad. - NoseCone.plots : plots.aero_surface_plots._NoseConePlots - This contains all the plots methods. Use help(NoseCone.plots) to know - more about it. - NoseCone.prints : prints.aero_surface_prints._NoseConePrints - This contains all the prints methods. Use help(NoseCone.prints) to know - more about it. - """ - - def __init__( # pylint: disable=too-many-statements - self, - length, - kind, - base_radius=None, - bluffness=None, - rocket_radius=None, - power=None, - name="Nose Cone", - ): - """Initializes the nose cone. It is used to define the nose cone - length, kind, center of pressure and lift coefficient curve. - - Parameters - ---------- - length : float - Nose cone length. Has units of length and must be given in meters. - kind : string - Nose cone kind. Can be "conical", "ogive", "elliptical", "tangent", - "von karman", "parabolic", "powerseries" or "lvhaack". If - "powerseries" is used, the "power" argument must be assigned to a - value between 0 and 1. - base_radius : float, optional - Nose cone base radius. Has units of length and must be given in - meters. If not given, the ratio between ``base_radius`` and - ``rocket_radius`` will be assumed as 1. - bluffness : float, optional - Ratio between the radius of the circle on the tip of the ogive and - the radius of the base of the ogive. Currently only used for the - nose cone's drawing. Must be between 0 and 1. Default is None, which - means that the nose cone will not have a sphere on the tip. If a - value is given, the nose cone's length will be reduced to account - for the addition of the sphere at the tip. Must be None or 0 if a - "powerseries" nose cone kind is specified. - rocket_radius : int, float, optional - The reference rocket radius used for lift coefficient normalization. - If not given, the ratio between ``base_radius`` and - ``rocket_radius`` will be assumed as 1. - power : float, optional - Factor that controls the bluntness of the shape for a power series - nose cone. Must be between 0 and 1. It is ignored when other nose - cone types are used. - name : str, optional - Nose cone name. Has no impact in simulation, as it is only used to - display data in a more organized matter. - - Returns - ------- - None - """ - super().__init__(name) - - self._rocket_radius = rocket_radius - self._base_radius = base_radius - self._length = length - if bluffness is not None: - if bluffness > 1 or bluffness < 0: - raise ValueError( - f"Bluffness ratio of {bluffness} is out of range. " - "It must be between 0 and 1." - ) - self._bluffness = bluffness - if kind == "powerseries": - # Checks if bluffness is not being used - if (self.bluffness is not None) and (self.bluffness != 0): - raise ValueError( - "Parameter 'bluffness' must be None or 0 when using a nose cone kind 'powerseries'." - ) - - if power is None: - raise ValueError( - "Parameter 'power' cannot be None when using a nose cone kind 'powerseries'." - ) - - if power > 1 or power <= 0: - raise ValueError( - f"Power value of {power} is out of range. It must be between 0 and 1." - ) - self._power = power - self.kind = kind - - self.evaluate_lift_coefficient() - self.evaluate_center_of_pressure() - - self.plots = _NoseConePlots(self) - self.prints = _NoseConePrints(self) - - @property - def rocket_radius(self): - return self._rocket_radius - - @rocket_radius.setter - def rocket_radius(self, value): - self._rocket_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_lift_coefficient() - self.evaluate_nose_shape() - - @property - def base_radius(self): - return self._base_radius - - @base_radius.setter - def base_radius(self, value): - self._base_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_lift_coefficient() - self.evaluate_nose_shape() - - @property - def length(self): - return self._length - - @length.setter - def length(self, value): - self._length = value - self.evaluate_center_of_pressure() - self.evaluate_nose_shape() - - @property - def power(self): - return self._power - - @power.setter - def power(self, value): - if value is not None: - if value > 1 or value <= 0: - raise ValueError( - f"Power value of {value} is out of range. It must be between 0 and 1." - ) - self._power = value - self.evaluate_k() - self.evaluate_center_of_pressure() - self.evaluate_nose_shape() - - @property - def kind(self): - return self._kind - - @kind.setter - def kind(self, value): # pylint: disable=too-many-statements - # Analyzes nosecone type - # Sets the k for Cp calculation - # Sets the function which creates the respective curve - self._kind = value - value = (value.replace(" ", "")).lower() - - if value == "conical": - self.k = 2 / 3 - self.y_nosecone = Function(lambda x: x * self.base_radius / self.length) - - elif value == "lvhaack": - self.k = 0.563 - - def theta(x): - return np.arccos(1 - 2 * max(min(x / self.length, 1), 0)) - - self.y_nosecone = Function( - lambda x: self.base_radius - * (theta(x) - np.sin(2 * theta(x)) / 2 + (np.sin(theta(x)) ** 3) / 3) - ** (0.5) - / (np.pi**0.5) - ) - - elif value in ["tangent", "tangentogive", "ogive"]: - rho = (self.base_radius**2 + self.length**2) / (2 * self.base_radius) - volume = np.pi * ( - self.length * rho**2 - - (self.length**3) / 3 - - (rho - self.base_radius) * rho**2 * np.arcsin(self.length / rho) - ) - area = np.pi * self.base_radius**2 - self.k = 1 - volume / (area * self.length) - self.y_nosecone = Function( - lambda x: np.sqrt(rho**2 - (min(x - self.length, 0)) ** 2) - + (self.base_radius - rho) - ) - - elif value == "elliptical": - self.k = 1 / 3 - self.y_nosecone = Function( - lambda x: self.base_radius - * np.sqrt(1 - ((x - self.length) / self.length) ** 2) - ) - - elif value == "vonkarman": - self.k = 0.5 - - def theta(x): - return np.arccos(1 - 2 * max(min(x / self.length, 1), 0)) - - self.y_nosecone = Function( - lambda x: self.base_radius - * (theta(x) - np.sin(2 * theta(x)) / 2) ** (0.5) - / (np.pi**0.5) - ) - elif value == "parabolic": - self.k = 0.5 - self.y_nosecone = Function( - lambda x: self.base_radius - * ((2 * x / self.length - (x / self.length) ** 2) / (2 - 1)) - ) - elif value == "powerseries": - self.k = (2 * self.power) / ((2 * self.power) + 1) - self.y_nosecone = Function( - lambda x: self.base_radius * np.power(x / self.length, self.power) - ) - else: - raise ValueError( - f"Nose Cone kind '{self.kind}' not found, " - + "please use one of the following Nose Cone kinds:" - + '\n\t"conical"' - + '\n\t"ogive"' - + '\n\t"lvhaack"' - + '\n\t"tangent"' - + '\n\t"vonkarman"' - + '\n\t"elliptical"' - + '\n\t"powerseries"' - + '\n\t"parabolic"\n' - ) - - self.evaluate_center_of_pressure() - self.evaluate_geometrical_parameters() - self.evaluate_nose_shape() - - @property - def bluffness(self): - return self._bluffness - - @bluffness.setter - def bluffness(self, value): - # prevents from setting bluffness on "powerseries" nose cones - if self.kind == "powerseries": - # Checks if bluffness is not being used - if (value is not None) and (value != 0): - raise ValueError( - "Parameter 'bluffness' must be None or 0 when using a nose cone kind 'powerseries'." - ) - if value is not None: - if value > 1 or value < 0: - raise ValueError( - f"Bluffness ratio of {value} is out of range. " - "It must be between 0 and 1." - ) - self._bluffness = value - self.evaluate_nose_shape() - - def evaluate_geometrical_parameters(self): - """Calculates and saves nose cone's radius ratio. - - Returns - ------- - None - """ - - # If base radius is not given, the ratio between base radius and - # rocket radius is assumed as 1, meaning that the nose cone has the - # same radius as the rocket - if self.base_radius is None and self.rocket_radius is not None: - self.radius_ratio = 1 - self.base_radius = self.rocket_radius - elif self.base_radius is not None and self.rocket_radius is None: - self.radius_ratio = 1 - self.rocket_radius = self.base_radius - # If base radius is given, the ratio between base radius and rocket - # radius is calculated - elif self.base_radius is not None and self.rocket_radius is not None: - self.radius_ratio = self.base_radius / self.rocket_radius - else: - raise ValueError( - "Either base radius or rocket radius must be given to " - "calculate the nose cone radius ratio." - ) - - self.fineness_ratio = self.length / (2 * self.base_radius) - - def evaluate_nose_shape(self): # pylint: disable=too-many-statements - """Calculates and saves nose cone's shape as lists and re-evaluates the - nose cone's length for a given bluffness ratio. The shape is saved as - two vectors, one for the x coordinates and one for the y coordinates. - - Returns - ------- - None - """ - number_of_points = 127 - density_modifier = 3 # increase density of points to improve accuracy - - def find_x_intercept(x): - # find the tangential intersection point between the circle and nosec curve - return x + self.y_nosecone(x) * self.y_nosecone.differentiate_complex_step( - x - ) - - # Calculate a function to find the radius of the nosecone curve - def find_radius(x): - return (self.y_nosecone(x) ** 2 + (x - find_x_intercept(x)) ** 2) ** 0.5 - - # Check bluffness circle and choose whether to use it or not - if self.bluffness is None or self.bluffness == 0: - # Set up parameters to continue without bluffness - r_circle, circle_center, x_init = 0, 0, 0 - else: - # Calculate circle radius - r_circle = self.bluffness * self.base_radius - if self.kind == "elliptical": - - def test_circle(x): - # set up a circle at the starting position to test bluffness - return np.sqrt(r_circle**2 - (x - r_circle) ** 2) - - # Check if bluffness circle is too small - if test_circle(1e-03) <= self.y_nosecone(1e-03): - # Raise a warning - warnings.warn( - "WARNING: The chosen bluffness ratio is too small for " - "the selected nose cone category, thereby the effective " - "bluffness will be 0." - ) - # Set up parameters to continue without bluffness - r_circle, circle_center, x_init = 0, 0, 0 - else: - # Find the intersection point between circle and nosecone curve - x_init = fsolve(lambda x: find_radius(x[0]) - r_circle, r_circle)[0] - circle_center = find_x_intercept(x_init) - else: - # Find the intersection point between circle and nosecone curve - x_init = fsolve(lambda x: find_radius(x[0]) - r_circle, r_circle)[0] - circle_center = find_x_intercept(x_init) - - # Calculate a function to create the circle at the correct position - def create_circle(x): - return abs(r_circle**2 - (x - circle_center) ** 2) ** 0.5 - - # Define a function for the final shape of the curve with a circle at the tip - def final_shape(x): - return self.y_nosecone(x) if x >= x_init else create_circle(x) - - # Vectorize the final_shape function - final_shape_vec = np.vectorize(final_shape) - - # Create the vectors X and Y with the points of the curve - nosecone_x = (self.length - (circle_center - r_circle)) * ( - np.linspace(0, 1, number_of_points) ** density_modifier - ) - nosecone_y = final_shape_vec(nosecone_x + (circle_center - r_circle)) - - # Evaluate final geometry parameters - self.shape_vec = [nosecone_x, nosecone_y] - if abs(nosecone_x[-1] - self.length) >= 0.001: # 1 millimeter - self._length = nosecone_x[-1] - print( - "Due to the chosen bluffness ratio, the nose " - f"cone length was reduced to {self.length} m." - ) - self.fineness_ratio = self.length / (2 * self.base_radius) - - def evaluate_lift_coefficient(self): - """Calculates and returns nose cone's lift coefficient. - The lift coefficient is saved and returned. This function - also calculates and saves its lift coefficient derivative. - - Returns - ------- - None - """ - # Calculate clalpha - # clalpha is currently a constant, meaning it is independent of Mach - # number. This is only valid for subsonic speeds. - # It must be set as a Function because it will be called and treated - # as a function of mach in the simulation. - self.clalpha = Function( - lambda mach: 2 / self._beta(mach) * self.radius_ratio**2, - "Mach", - f"Lift coefficient derivative for {self.name}", - ) - self.cl = Function( - lambda alpha, mach: self.clalpha(mach) * alpha, - ["Alpha (rad)", "Mach"], - "Cl", - ) - - def evaluate_k(self): - """Updates the self.k attribute used to compute the center of - pressure when using "powerseries" nose cones. - - Returns - ------- - None - """ - if self.kind == "powerseries": - self.k = (2 * self.power) / ((2 * self.power) + 1) - - def evaluate_center_of_pressure(self): - """Calculates and returns the center of pressure of the nose cone in - local coordinates. The center of pressure position is saved and stored - as a tuple. Local coordinate origin is found at the tip of the nose - cone. - - Returns - ------- - self.cp : tuple - Tuple containing cpx, cpy, cpz. - """ - - self.cpz = self.k * self.length - self.cpy = 0 - self.cpx = 0 - self.cp = (self.cpx, self.cpy, self.cpz) - return self.cp - - def draw(self): - """Draw the fin shape along with some important information, including - the center line, the quarter line and the center of pressure position. - - Returns - ------- - None - """ - self.plots.draw() - - def info(self): - """Prints and plots summarized information of the nose cone. - - Return - ------ - None - """ - self.prints.geometry() - self.prints.lift() - - def all_info(self): - """Prints and plots all the available information of the nose cone. - - Returns - ------- - None - """ - self.prints.all() - self.plots.all() - - -class Fins(AeroSurface): - """Abstract class that holds common methods for the fin classes. - Cannot be instantiated. - - Note - ---- - Local coordinate system: Z axis along the longitudinal axis of symmetry, - positive downwards (top -> bottom). Origin located at the top of the root - chord. - - Attributes - ---------- - Fins.n : int - Number of fins in fin set. - Fins.rocket_radius : float - The reference rocket radius used for lift coefficient normalization, - in meters. - Fins.airfoil : tuple - Tuple of two items. First is the airfoil lift curve. - Second is the unit of the curve (radians or degrees). - Fins.cant_angle : float - Fins cant angle with respect to the rocket centerline, in degrees. - Fins.changing_attribute_dict : dict - Dictionary that stores the name and the values of the attributes that - may be changed during a simulation. Useful for control systems. - Fins.cant_angle_rad : float - Fins cant angle with respect to the rocket centerline, in radians. - Fins.root_chord : float - Fin root chord in meters. - Fins.tip_chord : float - Fin tip chord in meters. - Fins.span : float - Fin span in meters. - Fins.name : string - Name of fin set. - Fins.sweep_length : float - Fins sweep length in meters. By sweep length, understand the axial - distance between the fin root leading edge and the fin tip leading edge - measured parallel to the rocket centerline. - Fins.sweep_angle : float - Fins sweep angle with respect to the rocket centerline. Must - be given in degrees. - Fins.d : float - Reference diameter of the rocket. Has units of length and is given - in meters. - Fins.ref_area : float - Reference area of the rocket. - Fins.Af : float - Area of the longitudinal section of each fin in the set. - Fins.AR : float - Aspect ratio of each fin in the set. - Fins.gamma_c : float - Fin mid-chord sweep angle. - Fins.Yma : float - Span wise position of the mean aerodynamic chord. - Fins.roll_geometrical_constant : float - Geometrical constant used in roll calculations. - Fins.tau : float - Geometrical relation used to simplify lift and roll calculations. - Fins.lift_interference_factor : float - Factor of Fin-Body interference in the lift coefficient. - Fins.cp : tuple - Tuple with the x, y and z local coordinates of the fin set center of - pressure. Has units of length and is given in meters. - Fins.cpx : float - Fin set local center of pressure x coordinate. Has units of length and - is given in meters. - Fins.cpy : float - Fin set local center of pressure y coordinate. Has units of length and - is given in meters. - Fins.cpz : float - Fin set local center of pressure z coordinate. Has units of length and - is given in meters. - Fins.cl : Function - Function which defines the lift coefficient as a function of the angle - of attack and the Mach number. Takes as input the angle of attack in - radians and the Mach number. Returns the lift coefficient. - Fins.clalpha : float - Lift coefficient slope. Has units of 1/rad. - Fins.roll_parameters : list - List containing the roll moment lift coefficient, the roll moment - damping coefficient and the cant angle in radians. - """ - - def __init__( - self, - n, - root_chord, - span, - rocket_radius, - cant_angle=0, - airfoil=None, - name="Fins", - ): - """Initialize Fins class. - - Parameters - ---------- - n : int - Number of fins, from 2 to infinity. - root_chord : int, float - Fin root chord in meters. - span : int, float - Fin span in meters. - rocket_radius : int, float - Reference rocket radius used for lift coefficient normalization. - cant_angle : int, float, optional - Fins cant angle with respect to the rocket centerline. Must - be given in degrees. - airfoil : tuple, optional - Default is null, in which case fins will be treated as flat plates. - Otherwise, if tuple, fins will be considered as airfoils. The - tuple's first item specifies the airfoil's lift coefficient - by angle of attack and must be either a .csv, .txt, ndarray - or callable. The .csv and .txt files can contain a single line - header and the first column must specify the angle of attack, while - the second column must specify the lift coefficient. The - ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] - where x0 is the angle of attack and y0 is the lift coefficient. - If callable, it should take an angle of attack as input and - return the lift coefficient at that angle of attack. - The tuple's second item is the unit of the angle of attack, - accepting either "radians" or "degrees". - name : str - Name of fin set. - - Returns - ------- - None - """ - - super().__init__(name) - - # Compute auxiliary geometrical parameters - d = 2 * rocket_radius - ref_area = np.pi * rocket_radius**2 # Reference area - - # Store values - self._n = n - self._rocket_radius = rocket_radius - self._airfoil = airfoil - self._cant_angle = cant_angle - self._root_chord = root_chord - self._span = span - self.name = name - self.d = d - self.ref_area = ref_area # Reference area - - @property - def n(self): - return self._n - - @n.setter - def n(self, value): - self._n = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def root_chord(self): - return self._root_chord - - @root_chord.setter - def root_chord(self, value): - self._root_chord = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def span(self): - return self._span - - @span.setter - def span(self, value): - self._span = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def rocket_radius(self): - return self._rocket_radius - - @rocket_radius.setter - def rocket_radius(self, value): - self._rocket_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def cant_angle(self): - return self._cant_angle - - @cant_angle.setter - def cant_angle(self, value): - self._cant_angle = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def airfoil(self): - return self._airfoil - - @airfoil.setter - def airfoil(self, value): - self._airfoil = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - def evaluate_lift_coefficient(self): - """Calculates and returns the fin set's lift coefficient. - The lift coefficient is saved and returned. This function - also calculates and saves the lift coefficient derivative - for a single fin and the lift coefficient derivative for - a number of n fins corrected for Fin-Body interference. - - Returns - ------- - None - """ - if not self.airfoil: - # Defines clalpha2D as 2*pi for planar fins - clalpha2D_incompressible = 2 * np.pi # pylint: disable=invalid-name - else: - # Defines clalpha2D as the derivative of the lift coefficient curve - # for the specific airfoil - self.airfoil_cl = Function( - self.airfoil[0], - interpolation="linear", - ) - - # Differentiating at alpha = 0 to get cl_alpha - clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step( # pylint: disable=invalid-name - x=1e-3, dx=1e-3 - ) - - # Convert to radians if needed - if self.airfoil[1] == "degrees": - clalpha2D_incompressible *= 180 / np.pi # pylint: disable=invalid-name - - # Correcting for compressible flow (apply Prandtl-Glauert correction) - clalpha2D = Function(lambda mach: clalpha2D_incompressible / self._beta(mach)) - - # Diederich's Planform Correlation Parameter - planform_correlation_parameter = ( - 2 * np.pi * self.AR / (clalpha2D * np.cos(self.gamma_c)) - ) # pylint: disable=invalid-name - - # Lift coefficient derivative for a single fin - def lift_source(mach): - return ( - clalpha2D(mach) - * planform_correlation_parameter(mach) - * (self.Af / self.ref_area) - * np.cos(self.gamma_c) - ) / ( - 2 - + planform_correlation_parameter(mach) - * np.sqrt(1 + (2 / planform_correlation_parameter(mach)) ** 2) - ) - - self.clalpha_single_fin = Function( - lift_source, - "Mach", - "Lift coefficient derivative for a single fin", - ) - - # Lift coefficient derivative for n fins corrected with Fin-Body interference - self.clalpha_multiple_fins = ( - self.lift_interference_factor - * self.fin_num_correction(self.n) - * self.clalpha_single_fin - ) # Function of mach number - self.clalpha_multiple_fins.set_inputs("Mach") - self.clalpha_multiple_fins.set_outputs( - f"Lift coefficient derivative for {self.n:.0f} fins" - ) - self.clalpha = self.clalpha_multiple_fins - - # Cl = clalpha * alpha - self.cl = Function( - lambda alpha, mach: alpha * self.clalpha_multiple_fins(mach), - ["Alpha (rad)", "Mach"], - "Lift coefficient", - ) - - return self.cl - - def evaluate_roll_parameters(self): - """Calculates and returns the fin set's roll coefficients. - The roll coefficients are saved in a list. - - Returns - ------- - self.roll_parameters : list - List containing the roll moment lift coefficient, the - roll moment damping coefficient and the cant angle in - radians - """ - - self.cant_angle_rad = np.radians(self.cant_angle) - - clf_delta = ( - self.roll_forcing_interference_factor - * self.n - * (self.Yma + self.rocket_radius) - * self.clalpha_single_fin - / self.d - ) # Function of mach number - clf_delta.set_inputs("Mach") - clf_delta.set_outputs("Roll moment forcing coefficient derivative") - cld_omega = ( - 2 - * self.roll_damping_interference_factor - * self.n - * self.clalpha_single_fin - * np.cos(self.cant_angle_rad) - * self.roll_geometrical_constant - / (self.ref_area * self.d**2) - ) # Function of mach number - cld_omega.set_inputs("Mach") - cld_omega.set_outputs("Roll moment damping coefficient derivative") - self.roll_parameters = [clf_delta, cld_omega, self.cant_angle_rad] - return self.roll_parameters - - @staticmethod - def fin_num_correction(n): - """Calculates a correction factor for the lift coefficient of multiple - fins. - The specifics values are documented at: - Niskanen, S. (2013). “OpenRocket technical documentation”. - In: Development of an Open Source model rocket simulation software. - - Parameters - ---------- - n : int - Number of fins. - - Returns - ------- - Corrector factor : int - Factor that accounts for the number of fins. - """ - corrector_factor = [2.37, 2.74, 2.99, 3.24] - if 5 <= n <= 8: - return corrector_factor[n - 5] - else: - return n / 2 - - def draw(self): - """Draw the fin shape along with some important information, including - the center line, the quarter line and the center of pressure position. - - Returns - ------- - None - """ - self.plots.draw() - - -class TrapezoidalFins(Fins): - """Class that defines and holds information for a trapezoidal fin set. - - This class inherits from the Fins class. - - Note - ---- - Local coordinate system: Z axis along the longitudinal axis of symmetry, - positive downwards (top -> bottom). Origin located at the top of the root - chord. - - See Also - -------- - Fins - - Attributes - ---------- - TrapezoidalFins.n : int - Number of fins in fin set. - TrapezoidalFins.rocket_radius : float - The reference rocket radius used for lift coefficient normalization, in - meters. - TrapezoidalFins.airfoil : tuple - Tuple of two items. First is the airfoil lift curve. - Second is the unit of the curve (radians or degrees). - TrapezoidalFins.cant_angle : float - Fins cant angle with respect to the rocket centerline, in degrees. - TrapezoidalFins.changing_attribute_dict : dict - Dictionary that stores the name and the values of the attributes that - may be changed during a simulation. Useful for control systems. - TrapezoidalFins.cant_angle_rad : float - Fins cant angle with respect to the rocket centerline, in radians. - TrapezoidalFins.root_chord : float - Fin root chord in meters. - TrapezoidalFins.tip_chord : float - Fin tip chord in meters. - TrapezoidalFins.span : float - Fin span in meters. - TrapezoidalFins.name : string - Name of fin set. - TrapezoidalFins.sweep_length : float - Fins sweep length in meters. By sweep length, understand the axial - distance between the fin root leading edge and the fin tip leading edge - measured parallel to the rocket centerline. - TrapezoidalFins.sweep_angle : float - Fins sweep angle with respect to the rocket centerline. Must - be given in degrees. - TrapezoidalFins.d : float - Reference diameter of the rocket, in meters. - TrapezoidalFins.ref_area : float - Reference area of the rocket, in m². - TrapezoidalFins.Af : float - Area of the longitudinal section of each fin in the set. - TrapezoidalFins.AR : float - Aspect ratio of each fin in the set - TrapezoidalFins.gamma_c : float - Fin mid-chord sweep angle. - TrapezoidalFins.Yma : float - Span wise position of the mean aerodynamic chord. - TrapezoidalFins.roll_geometrical_constant : float - Geometrical constant used in roll calculations. - TrapezoidalFins.tau : float - Geometrical relation used to simplify lift and roll calculations. - TrapezoidalFins.lift_interference_factor : float - Factor of Fin-Body interference in the lift coefficient. - TrapezoidalFins.cp : tuple - Tuple with the x, y and z local coordinates of the fin set center of - pressure. Has units of length and is given in meters. - TrapezoidalFins.cpx : float - Fin set local center of pressure x coordinate. Has units of length and - is given in meters. - TrapezoidalFins.cpy : float - Fin set local center of pressure y coordinate. Has units of length and - is given in meters. - TrapezoidalFins.cpz : float - Fin set local center of pressure z coordinate. Has units of length and - is given in meters. - TrapezoidalFins.cl : Function - Function which defines the lift coefficient as a function of the angle - of attack and the Mach number. Takes as input the angle of attack in - radians and the Mach number. Returns the lift coefficient. - TrapezoidalFins.clalpha : float - Lift coefficient slope. Has units of 1/rad. - """ - - def __init__( - self, - n, - root_chord, - tip_chord, - span, - rocket_radius, - cant_angle=0, - sweep_length=None, - sweep_angle=None, - airfoil=None, - name="Fins", - ): - """Initialize TrapezoidalFins class. - - Parameters - ---------- - n : int - Number of fins, from 2 to infinity. - root_chord : int, float - Fin root chord in meters. - tip_chord : int, float - Fin tip chord in meters. - span : int, float - Fin span in meters. - rocket_radius : int, float - Reference radius to calculate lift coefficient, in meters. - cant_angle : int, float, optional - Fins cant angle with respect to the rocket centerline. Must - be given in degrees. - sweep_length : int, float, optional - Fins sweep length in meters. By sweep length, understand the axial - distance between the fin root leading edge and the fin tip leading - edge measured parallel to the rocket centerline. If not given, the - sweep length is assumed to be equal the root chord minus the tip - chord, in which case the fin is a right trapezoid with its base - perpendicular to the rocket's axis. Cannot be used in conjunction - with sweep_angle. - sweep_angle : int, float, optional - Fins sweep angle with respect to the rocket centerline. Must - be given in degrees. If not given, the sweep angle is automatically - calculated, in which case the fin is assumed to be a right trapezoid - with its base perpendicular to the rocket's axis. - Cannot be used in conjunction with sweep_length. - airfoil : tuple, optional - Default is null, in which case fins will be treated as flat plates. - Otherwise, if tuple, fins will be considered as airfoils. The - tuple's first item specifies the airfoil's lift coefficient - by angle of attack and must be either a .csv, .txt, ndarray - or callable. The .csv and .txt files can contain a single line - header and the first column must specify the angle of attack, while - the second column must specify the lift coefficient. The - ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] - where x0 is the angle of attack and y0 is the lift coefficient. - If callable, it should take an angle of attack as input and - return the lift coefficient at that angle of attack. - The tuple's second item is the unit of the angle of attack, - accepting either "radians" or "degrees". - name : str - Name of fin set. - - Returns - ------- - None - """ - - super().__init__( - n, - root_chord, - span, - rocket_radius, - cant_angle, - airfoil, - name, - ) - - # Check if sweep angle or sweep length is given - if sweep_length is not None and sweep_angle is not None: - raise ValueError("Cannot use sweep_length and sweep_angle together") - elif sweep_angle is not None: - sweep_length = np.tan(sweep_angle * np.pi / 180) * span - elif sweep_length is None: - sweep_length = root_chord - tip_chord - else: - # Sweep length is given - pass - - self._tip_chord = tip_chord - self._sweep_length = sweep_length - self._sweep_angle = sweep_angle - - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - self.prints = _TrapezoidalFinsPrints(self) - self.plots = _TrapezoidalFinsPlots(self) - - @property - def tip_chord(self): - return self._tip_chord - - @tip_chord.setter - def tip_chord(self, value): - self._tip_chord = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def sweep_angle(self): - return self._sweep_angle - - @sweep_angle.setter - def sweep_angle(self, value): - self._sweep_angle = value - self._sweep_length = np.tan(value * np.pi / 180) * self.span - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def sweep_length(self): - return self._sweep_length - - @sweep_length.setter - def sweep_length(self, value): - self._sweep_length = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - def evaluate_center_of_pressure(self): - """Calculates and returns the center of pressure of the fin set in local - coordinates. The center of pressure position is saved and stored as a - tuple. - - Returns - ------- - None - """ - # Center of pressure position in local coordinates - cpz = (self.sweep_length / 3) * ( - (self.root_chord + 2 * self.tip_chord) / (self.root_chord + self.tip_chord) - ) + (1 / 6) * ( - self.root_chord - + self.tip_chord - - self.root_chord * self.tip_chord / (self.root_chord + self.tip_chord) - ) - self.cpx = 0 - self.cpy = 0 - self.cpz = cpz - self.cp = (self.cpx, self.cpy, self.cpz) - - def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements - """Calculates and saves fin set's geometrical parameters such as the - fins' area, aspect ratio and parameters for roll movement. - - Returns - ------- - None - """ - # pylint: disable=invalid-name - Yr = self.root_chord + self.tip_chord - Af = Yr * self.span / 2 # Fin area - AR = 2 * self.span**2 / Af # Fin aspect ratio - gamma_c = np.arctan( - (self.sweep_length + 0.5 * self.tip_chord - 0.5 * self.root_chord) - / (self.span) - ) - Yma = ( - (self.span / 3) * (self.root_chord + 2 * self.tip_chord) / Yr - ) # Span wise coord of mean aero chord - - # Fin–body interference correction parameters - tau = (self.span + self.rocket_radius) / self.rocket_radius - lift_interference_factor = 1 + 1 / tau - lambda_ = self.tip_chord / self.root_chord - - # Parameters for Roll Moment. - # Documented at: https://docs.rocketpy.org/en/latest/technical/ - roll_geometrical_constant = ( - (self.root_chord + 3 * self.tip_chord) * self.span**3 - + 4 - * (self.root_chord + 2 * self.tip_chord) - * self.rocket_radius - * self.span**2 - + 6 * (self.root_chord + self.tip_chord) * self.span * self.rocket_radius**2 - ) / 12 - roll_damping_interference_factor = 1 + ( - ((tau - lambda_) / (tau)) - ((1 - lambda_) / (tau - 1)) * np.log(tau) - ) / ( - ((tau + 1) * (tau - lambda_)) / (2) - - ((1 - lambda_) * (tau**3 - 1)) / (3 * (tau - 1)) - ) - roll_forcing_interference_factor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Store values - self.Yr = Yr - self.Af = Af # Fin area - self.AR = AR # Aspect Ratio - self.gamma_c = gamma_c # Mid chord angle - self.Yma = Yma # Span wise coord of mean aero chord - self.roll_geometrical_constant = roll_geometrical_constant - self.tau = tau - self.lift_interference_factor = lift_interference_factor - self.λ = lambda_ # pylint: disable=non-ascii-name - self.roll_damping_interference_factor = roll_damping_interference_factor - self.roll_forcing_interference_factor = roll_forcing_interference_factor - - self.evaluate_shape() - - def evaluate_shape(self): - if self.sweep_length: - points = [ - (0, 0), - (self.sweep_length, self.span), - (self.sweep_length + self.tip_chord, self.span), - (self.root_chord, 0), - ] - else: - points = [ - (0, 0), - (self.root_chord - self.tip_chord, self.span), - (self.root_chord, self.span), - (self.root_chord, 0), - ] - - x_array, y_array = zip(*points) - self.shape_vec = [np.array(x_array), np.array(y_array)] - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() - - -class EllipticalFins(Fins): - """Class that defines and holds information for an elliptical fin set. - - This class inherits from the Fins class. - - Note - ---- - Local coordinate system: Z axis along the longitudinal axis of symmetry, - positive downwards (top -> bottom). Origin located at the top of the root - chord. - - See Also - -------- - Fins - - Attributes - ---------- - EllipticalFins.n : int - Number of fins in fin set. - EllipticalFins.rocket_radius : float - The reference rocket radius used for lift coefficient normalization, in - meters. - EllipticalFins.airfoil : tuple - Tuple of two items. First is the airfoil lift curve. - Second is the unit of the curve (radians or degrees) - EllipticalFins.cant_angle : float - Fins cant angle with respect to the rocket centerline, in degrees. - EllipticalFins.changing_attribute_dict : dict - Dictionary that stores the name and the values of the attributes that - may be changed during a simulation. Useful for control systems. - EllipticalFins.cant_angle_rad : float - Fins cant angle with respect to the rocket centerline, in radians. - EllipticalFins.root_chord : float - Fin root chord in meters. - EllipticalFins.span : float - Fin span in meters. - EllipticalFins.name : string - Name of fin set. - EllipticalFins.sweep_length : float - Fins sweep length in meters. By sweep length, understand the axial - distance between the fin root leading edge and the fin tip leading edge - measured parallel to the rocket centerline. - EllipticalFins.sweep_angle : float - Fins sweep angle with respect to the rocket centerline. Must - be given in degrees. - EllipticalFins.d : float - Reference diameter of the rocket, in meters. - EllipticalFins.ref_area : float - Reference area of the rocket. - EllipticalFins.Af : float - Area of the longitudinal section of each fin in the set. - EllipticalFins.AR : float - Aspect ratio of each fin in the set. - EllipticalFins.gamma_c : float - Fin mid-chord sweep angle. - EllipticalFins.Yma : float - Span wise position of the mean aerodynamic chord. - EllipticalFins.roll_geometrical_constant : float - Geometrical constant used in roll calculations. - EllipticalFins.tau : float - Geometrical relation used to simplify lift and roll calculations. - EllipticalFins.lift_interference_factor : float - Factor of Fin-Body interference in the lift coefficient. - EllipticalFins.cp : tuple - Tuple with the x, y and z local coordinates of the fin set center of - pressure. Has units of length and is given in meters. - EllipticalFins.cpx : float - Fin set local center of pressure x coordinate. Has units of length and - is given in meters. - EllipticalFins.cpy : float - Fin set local center of pressure y coordinate. Has units of length and - is given in meters. - EllipticalFins.cpz : float - Fin set local center of pressure z coordinate. Has units of length and - is given in meters. - EllipticalFins.cl : Function - Function which defines the lift coefficient as a function of the angle - of attack and the Mach number. Takes as input the angle of attack in - radians and the Mach number. Returns the lift coefficient. - EllipticalFins.clalpha : float - Lift coefficient slope. Has units of 1/rad. - """ - - def __init__( - self, - n, - root_chord, - span, - rocket_radius, - cant_angle=0, - airfoil=None, - name="Fins", - ): - """Initialize EllipticalFins class. - - Parameters - ---------- - n : int - Number of fins, from 2 to infinity. - root_chord : int, float - Fin root chord in meters. - span : int, float - Fin span in meters. - rocket_radius : int, float - Reference radius to calculate lift coefficient, in meters. - cant_angle : int, float, optional - Fins cant angle with respect to the rocket centerline. Must - be given in degrees. - sweep_length : int, float, optional - Fins sweep length in meters. By sweep length, understand the axial - distance between the fin root leading edge and the fin tip leading - edge measured parallel to the rocket centerline. If not given, the - sweep length is assumed to be equal the root chord minus the tip - chord, in which case the fin is a right trapezoid with its base - perpendicular to the rocket's axis. Cannot be used in conjunction - with sweep_angle. - sweep_angle : int, float, optional - Fins sweep angle with respect to the rocket centerline. Must - be given in degrees. If not given, the sweep angle is automatically - calculated, in which case the fin is assumed to be a right trapezoid - with its base perpendicular to the rocket's axis. - Cannot be used in conjunction with sweep_length. - airfoil : tuple, optional - Default is null, in which case fins will be treated as flat plates. - Otherwise, if tuple, fins will be considered as airfoils. The - tuple's first item specifies the airfoil's lift coefficient - by angle of attack and must be either a .csv, .txt, ndarray - or callable. The .csv and .txt files can contain a single line - header and the first column must specify the angle of attack, while - the second column must specify the lift coefficient. The - ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] - where x0 is the angle of attack and y0 is the lift coefficient. - If callable, it should take an angle of attack as input and - return the lift coefficient at that angle of attack. - The tuple's second item is the unit of the angle of attack, - accepting either "radians" or "degrees". - name : str - Name of fin set. - - Returns - ------- - None - """ - - super().__init__( - n, - root_chord, - span, - rocket_radius, - cant_angle, - airfoil, - name, - ) - - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - self.prints = _EllipticalFinsPrints(self) - self.plots = _EllipticalFinsPlots(self) - - def evaluate_center_of_pressure(self): - """Calculates and returns the center of pressure of the fin set in local - coordinates. The center of pressure position is saved and stored as a - tuple. - - Returns - ------- - None - """ - # Center of pressure position in local coordinates - cpz = 0.288 * self.root_chord - self.cpx = 0 - self.cpy = 0 - self.cpz = cpz - self.cp = (self.cpx, self.cpy, self.cpz) - - def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements - """Calculates and saves fin set's geometrical parameters such as the - fins' area, aspect ratio and parameters for roll movement. - - Returns - ------- - None - """ - - # Compute auxiliary geometrical parameters - Af = ( # Fin area # pylint: disable=invalid-name - np.pi * self.root_chord / 2 * self.span - ) / 2 - gamma_c = 0 # Zero for elliptical fins - AR = 2 * self.span**2 / Af # Fin aspect ratio # pylint: disable=invalid-name - Yma = ( # pylint: disable=invalid-name - self.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) - ) # Span wise coord of mean aero chord - roll_geometrical_constant = ( - self.root_chord - * self.span - * ( - 3 * np.pi * self.span**2 - + 32 * self.rocket_radius * self.span - + 12 * np.pi * self.rocket_radius**2 - ) - / 48 - ) - - # Fin–body interference correction parameters - tau = (self.span + self.rocket_radius) / self.rocket_radius - lift_interference_factor = 1 + 1 / tau - if self.span > self.rocket_radius: - roll_damping_interference_factor = 1 + ( - (self.rocket_radius**2) - * ( - 2 - * (self.rocket_radius**2) - * np.sqrt(self.span**2 - self.rocket_radius**2) - * np.log( - ( - 2 - * self.span - * np.sqrt(self.span**2 - self.rocket_radius**2) - + 2 * self.span**2 - ) - / self.rocket_radius - ) - - 2 - * (self.rocket_radius**2) - * np.sqrt(self.span**2 - self.rocket_radius**2) - * np.log(2 * self.span) - + 2 * self.span**3 - - np.pi * self.rocket_radius * self.span**2 - - 2 * (self.rocket_radius**2) * self.span - + np.pi * self.rocket_radius**3 - ) - ) / ( - 2 - * (self.span**2) - * (self.span / 3 + np.pi * self.rocket_radius / 4) - * (self.span**2 - self.rocket_radius**2) - ) - elif self.span < self.rocket_radius: - roll_damping_interference_factor = 1 - ( - self.rocket_radius**2 - * ( - 2 * self.span**3 - - np.pi * self.span**2 * self.rocket_radius - - 2 * self.span * self.rocket_radius**2 - + np.pi * self.rocket_radius**3 - + 2 - * self.rocket_radius**2 - * np.sqrt(-self.span**2 + self.rocket_radius**2) - * np.arctan( - (self.span) / (np.sqrt(-self.span**2 + self.rocket_radius**2)) - ) - - np.pi - * self.rocket_radius**2 - * np.sqrt(-self.span**2 + self.rocket_radius**2) - ) - ) / ( - 2 - * self.span - * (-self.span**2 + self.rocket_radius**2) - * (self.span**2 / 3 + np.pi * self.span * self.rocket_radius / 4) - ) - else: - roll_damping_interference_factor = (28 - 3 * np.pi) / (4 + 3 * np.pi) - - roll_forcing_interference_factor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Store values - self.Af = Af # Fin area # pylint: disable=invalid-name - self.AR = AR # Fin aspect ratio # pylint: disable=invalid-name - self.gamma_c = gamma_c # Mid chord angle - self.Yma = ( # pylint: disable=invalid-name - Yma # Span wise coord of mean aero chord # pylint: disable=invalid-name - ) - self.roll_geometrical_constant = roll_geometrical_constant - self.tau = tau - self.lift_interference_factor = lift_interference_factor - self.roll_damping_interference_factor = roll_damping_interference_factor - self.roll_forcing_interference_factor = roll_forcing_interference_factor - - self.evaluate_shape() - - def evaluate_shape(self): - angles = np.arange(0, 180, 5) - x_array = self.root_chord / 2 + self.root_chord / 2 * np.cos(np.radians(angles)) - y_array = self.span * np.sin(np.radians(angles)) - self.shape_vec = [x_array, y_array] - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() - - -class Tail(AeroSurface): - """Class that defines a tail. Currently only accepts conical tails. - - Note - ---- - Local coordinate system: Z axis along the longitudinal axis of symmetry, - positive downwards (top -> bottom). Origin located at top of the tail - (generally the portion closest to the rocket's nose). - - Attributes - ---------- - Tail.top_radius : int, float - Radius of the top of the tail. The top radius is defined as the radius - of the transversal section that is closest to the rocket's nose. - Tail.bottom_radius : int, float - Radius of the bottom of the tail. - Tail.length : int, float - Length of the tail. The length is defined as the distance between the - top and bottom of the tail. The length is measured along the rocket's - longitudinal axis. Has the unit of meters. - Tail.rocket_radius: int, float - The reference rocket radius used for lift coefficient normalization in - meters. - Tail.name : str - Name of the tail. Default is 'Tail'. - Tail.cpx : int, float - x local coordinate of the center of pressure of the tail. - Tail.cpy : int, float - y local coordinate of the center of pressure of the tail. - Tail.cpz : int, float - z local coordinate of the center of pressure of the tail. - Tail.cp : tuple - Tuple containing the coordinates of the center of pressure of the tail. - Tail.cl : Function - Function that returns the lift coefficient of the tail. The function - is defined as a function of the angle of attack and the mach number. - Tail.clalpha : float - Lift coefficient slope. Has the unit of 1/rad. - Tail.slant_length : float - Slant length of the tail. The slant length is defined as the distance - between the top and bottom of the tail. The slant length is measured - along the tail's slant axis. Has the unit of meters. - Tail.surface_area : float - Surface area of the tail. Has the unit of meters squared. - """ - - def __init__(self, top_radius, bottom_radius, length, rocket_radius, name="Tail"): - """Initializes the tail object by computing and storing the most - important values. - - Parameters - ---------- - top_radius : int, float - Radius of the top of the tail. The top radius is defined as the - radius of the transversal section that is closest to the rocket's - nose. - bottom_radius : int, float - Radius of the bottom of the tail. - length : int, float - Length of the tail. - rocket_radius : int, float - The reference rocket radius used for lift coefficient normalization. - name : str - Name of the tail. Default is 'Tail'. - - Returns - ------- - None - """ - super().__init__(name) - - self._top_radius = top_radius - self._bottom_radius = bottom_radius - self._length = length - self._rocket_radius = rocket_radius - - self.evaluate_geometrical_parameters() - self.evaluate_lift_coefficient() - self.evaluate_center_of_pressure() - - self.plots = _TailPlots(self) - self.prints = _TailPrints(self) - - @property - def top_radius(self): - return self._top_radius - - @top_radius.setter - def top_radius(self, value): - self._top_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_lift_coefficient() - self.evaluate_center_of_pressure() - - @property - def bottom_radius(self): - return self._bottom_radius - - @bottom_radius.setter - def bottom_radius(self, value): - self._bottom_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_lift_coefficient() - self.evaluate_center_of_pressure() - - @property - def length(self): - return self._length - - @length.setter - def length(self, value): - self._length = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - - @property - def rocket_radius(self): - return self._rocket_radius - - @rocket_radius.setter - def rocket_radius(self, value): - self._rocket_radius = value - self.evaluate_lift_coefficient() - - def evaluate_geometrical_parameters(self): - """Calculates and saves tail's slant length and surface area. - - Returns - ------- - None - """ - self.slant_length = np.sqrt( - (self.length) ** 2 + (self.top_radius - self.bottom_radius) ** 2 - ) - self.surface_area = ( - np.pi * self.slant_length * (self.top_radius + self.bottom_radius) - ) - self.evaluate_shape() - - def evaluate_shape(self): - # Assuming the tail is a cone, calculate the shape vector - self.shape_vec = [ - np.array([0, self.length]), - np.array([self.top_radius, self.bottom_radius]), - ] - - def evaluate_lift_coefficient(self): - """Calculates and returns tail's lift coefficient. - The lift coefficient is saved and returned. This function - also calculates and saves its lift coefficient derivative. - - Returns - ------- - None - """ - # Calculate clalpha - # clalpha is currently a constant, meaning it is independent of Mach - # number. This is only valid for subsonic speeds. - # It must be set as a Function because it will be called and treated - # as a function of mach in the simulation. - self.clalpha = Function( - lambda mach: 2 - / self._beta(mach) - * ( - (self.bottom_radius / self.rocket_radius) ** 2 - - (self.top_radius / self.rocket_radius) ** 2 - ), - "Mach", - f"Lift coefficient derivative for {self.name}", - ) - self.cl = Function( - lambda alpha, mach: self.clalpha(mach) * alpha, - ["Alpha (rad)", "Mach"], - "Cl", - ) - - def evaluate_center_of_pressure(self): - """Calculates and returns the center of pressure of the tail in local - coordinates. The center of pressure position is saved and stored as a - tuple. - - Returns - ------- - None - """ - # Calculate cp position in local coordinates - r = self.top_radius / self.bottom_radius - cpz = (self.length / 3) * (1 + (1 - r) / (1 - r**2)) - - # Store values as class attributes - self.cpx = 0 - self.cpy = 0 - self.cpz = cpz - self.cp = (self.cpx, self.cpy, self.cpz) - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() - - -class RailButtons(AeroSurface): - """Class that defines a rail button pair or group. - - Attributes - ---------- - RailButtons.buttons_distance : int, float - Distance between the two rail buttons closest to the nozzle. - RailButtons.angular_position : int, float - Angular position of the rail buttons in degrees measured - as the rotation around the symmetry axis of the rocket - relative to one of the other principal axis. - """ - - def __init__(self, buttons_distance, angular_position=45, name="Rail Buttons"): - """Initializes RailButtons Class. - - Parameters - ---------- - buttons_distance : int, float - Distance between the first and the last rail button in meters. - angular_position : int, float, optional - Angular position of the rail buttons in degrees measured - as the rotation around the symmetry axis of the rocket - relative to one of the other principal axis. - name : string, optional - Name of the rail buttons. Default is "Rail Buttons". - - Returns - ------- - None - - """ - super().__init__(name) - self.buttons_distance = buttons_distance - self.angular_position = angular_position - self.name = name - - self.evaluate_lift_coefficient() - self.evaluate_center_of_pressure() - - self.prints = _RailButtonsPrints(self) - - def evaluate_center_of_pressure(self): - """Evaluates the center of pressure of the rail buttons. Rail buttons - do not contribute to the center of pressure of the rocket. - - Returns - ------- - None - """ - self.cpx = 0 - self.cpy = 0 - self.cpz = 0 - self.cp = (self.cpx, self.cpy, self.cpz) - - def evaluate_lift_coefficient(self): - """Evaluates the lift coefficient curve of the rail buttons. Rail - buttons do not contribute to the lift coefficient of the rocket. - - Returns - ------- - None - """ - self.clalpha = Function( - lambda mach: 0, - "Mach", - f"Lift coefficient derivative for {self.name}", - ) - self.cl = Function( - lambda alpha, mach: 0, - ["Alpha (rad)", "Mach"], - "Cl", - ) - - def evaluate_geometrical_parameters(self): - """Evaluates the geometrical parameters of the rail buttons. Rail - buttons do not contribute to the geometrical parameters of the rocket. - - Returns - ------- - None - """ - - def info(self): - """Prints out all the information about the Rail Buttons. - - Returns - ------- - None - """ - self.prints.geometry() - - def all_info(self): - """Returns all info of the Rail Buttons. - - Returns - ------- - None - """ - self.prints.all() - - -class AirBrakes(AeroSurface): - """AirBrakes class. Inherits from AeroSurface. - - Attributes - ---------- - AirBrakes.drag_coefficient : Function - Drag coefficient as a function of deployment level and Mach number. - AirBrakes.drag_coefficient_curve : int, float, callable, array, string, Function - Curve that defines the drag coefficient as a function of deployment level - and Mach number. Used as the source of `AirBrakes.drag_coefficient`. - AirBrakes.deployment_level : float - Current deployment level, ranging from 0 to 1. Deployment level is the - fraction of the total airbrake area that is deployed. - AirBrakes.reference_area : int, float - Reference area used to calculate the drag force of the air brakes - from the drag coefficient curve. Units of m^2. - AirBrakes.clamp : bool, optional - If True, the simulation will clamp the deployment level to 0 or 1 if - the deployment level is out of bounds. If False, the simulation will - not clamp the deployment level and will instead raise a warning if - the deployment level is out of bounds. Default is True. - AirBrakes.name : str - Name of the air brakes. - """ - - def __init__( - self, - drag_coefficient_curve, - reference_area, - clamp=True, - override_rocket_drag=False, - deployment_level=0, - name="AirBrakes", - ): - """Initializes the AirBrakes class. - - Parameters - ---------- - drag_coefficient_curve : int, float, callable, array, string, Function - This parameter represents the drag coefficient associated with the - air brakes and/or the entire rocket, depending on the value of - ``override_rocket_drag``. - - - If a constant, it should be an integer or a float representing a - fixed drag coefficient value. - - If a function, it must take two parameters: deployment level and - Mach number, and return the drag coefficient. This function allows - for dynamic computation based on deployment and Mach number. - - If an array, it should be a 2D array with three columns: the first - column for deployment level, the second for Mach number, and the - third for the corresponding drag coefficient. - - If a string, it should be the path to a .csv or .txt file. The - file must contain three columns: the first for deployment level, - the second for Mach number, and the third for the drag - coefficient. - - If a Function, it must take two parameters: deployment level and - Mach number, and return the drag coefficient. - - .. note:: For ``override_rocket_drag = False``, at - deployment level 0, the drag coefficient is assumed to be 0, - independent of the input drag coefficient curve. This means that - the simulation always considers that at a deployment level of 0, - the air brakes are completely retracted and do not contribute to - the drag of the rocket. - - reference_area : int, float - Reference area used to calculate the drag force of the air brakes - from the drag coefficient curve. Units of m^2. - clamp : bool, optional - If True, the simulation will clamp the deployment level to 0 or 1 if - the deployment level is out of bounds. If False, the simulation will - not clamp the deployment level and will instead raise a warning if - the deployment level is out of bounds. Default is True. - override_rocket_drag : bool, optional - If False, the air brakes drag coefficient will be added to the - rocket's power off drag coefficient curve. If True, during the - simulation, the rocket's power off drag will be ignored and the air - brakes drag coefficient will be used for the entire rocket instead. - Default is False. - deployment_level : float, optional - Initial deployment level, ranging from 0 to 1. Deployment level is - the fraction of the total airbrake area that is Deployment. Default - is 0. - name : str, optional - Name of the air brakes. Default is "AirBrakes". - - Returns - ------- - None - """ - super().__init__(name) - self.drag_coefficient_curve = drag_coefficient_curve - self.drag_coefficient = Function( - drag_coefficient_curve, - inputs=["Deployment Level", "Mach"], - outputs="Drag Coefficient", - ) - self.reference_area = reference_area - self.clamp = clamp - self.override_rocket_drag = override_rocket_drag - self.initial_deployment_level = deployment_level - self.deployment_level = deployment_level - self.prints = _AirBrakesPrints(self) - self.plots = _AirBrakesPlots(self) - - @property - def deployment_level(self): - """Returns the deployment level of the air brakes.""" - return self._deployment_level - - @deployment_level.setter - def deployment_level(self, value): - # Check if deployment level is within bounds and warn user if not - if value < 0 or value > 1: - # Clamp deployment level if clamp is True - if self.clamp: - # Make sure deployment level is between 0 and 1 - value = np.clip(value, 0, 1) - else: - # Raise warning if clamp is False - warnings.warn( - f"Deployment level of {self.name} is smaller than 0 or " - + "larger than 1. Extrapolation for the drag coefficient " - + "curve will be used." - ) - self._deployment_level = value - - def _reset(self): - """Resets the air brakes to their initial state. This is ran at the - beginning of each simulation to ensure the air brakes are in the correct - state.""" - self.deployment_level = self.initial_deployment_level - - def evaluate_center_of_pressure(self): - """Evaluates the center of pressure of the aerodynamic surface in local - coordinates. - - For air brakes, all components of the center of pressure position are - 0. - - Returns - ------- - None - """ - self.cpx = 0 - self.cpy = 0 - self.cpz = 0 - self.cp = (self.cpx, self.cpy, self.cpz) - - def evaluate_lift_coefficient(self): - """Evaluates the lift coefficient curve of the aerodynamic surface. - - For air brakes, the current model assumes no lift is generated. - Therefore, the lift coefficient (C_L) and its derivative relative to the - angle of attack (C_L_alpha), is 0. - - Returns - ------- - None - """ - self.clalpha = Function( - lambda mach: 0, - "Mach", - f"Lift coefficient derivative for {self.name}", - ) - self.cl = Function( - lambda alpha, mach: 0, - ["Alpha (rad)", "Mach"], - "Lift Coefficient", - ) - - def evaluate_geometrical_parameters(self): - """Evaluates the geometrical parameters of the aerodynamic surface. - - Returns - ------- - None - """ - - def info(self): - """Prints and plots summarized information of the aerodynamic surface. - - Returns - ------- - None - """ - self.prints.geometry() - - def all_info(self): - """Prints and plots all information of the aerodynamic surface. - - Returns - ------- - None - """ - self.info() - self.plots.drag_coefficient_curve() diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py new file mode 100644 index 000000000..9d9c68586 --- /dev/null +++ b/rocketpy/rocket/aero_surface/__init__.py @@ -0,0 +1,6 @@ +from rocketpy.rocket.aero_surface.aero_surface import AeroSurface +from rocketpy.rocket.aero_surface.air_brakes import AirBrakes +from rocketpy.rocket.aero_surface.fins import EllipticalFins, Fins, TrapezoidalFins +from rocketpy.rocket.aero_surface.nose_cone import NoseCone +from rocketpy.rocket.aero_surface.rail_buttons import RailButtons +from rocketpy.rocket.aero_surface.tail import Tail diff --git a/rocketpy/rocket/aero_surface/aero_surface.py b/rocketpy/rocket/aero_surface/aero_surface.py new file mode 100644 index 000000000..81b5610e3 --- /dev/null +++ b/rocketpy/rocket/aero_surface/aero_surface.py @@ -0,0 +1,92 @@ +from abc import ABC, abstractmethod + +import numpy as np + + +class AeroSurface(ABC): + """Abstract class used to define aerodynamic surfaces.""" + + def __init__(self, name, reference_area, reference_length): + self.reference_area = reference_area + self.reference_length = reference_length + self.name = name + + self.cpx = 0 + self.cpy = 0 + self.cpz = 0 + + @staticmethod + def _beta(mach): + """Defines a parameter that is often used in aerodynamic + equations. It is commonly used in the Prandtl factor which + corrects subsonic force coefficients for compressible flow. + This is applied to the lift coefficient of the nose cone, + fins and tails/transitions as in [1]. + + Parameters + ---------- + mach : int, float + Number of mach. + + Returns + ------- + beta : int, float + Value that characterizes flow speed based on the mach number. + + References + ---------- + [1] Barrowman, James S. https://arc.aiaa.org/doi/10.2514/6.1979-504 + """ + + if mach < 0.8: + return np.sqrt(1 - mach**2) + elif mach < 1.1: + return np.sqrt(1 - 0.8**2) + else: + return np.sqrt(mach**2 - 1) + + @abstractmethod + def evaluate_center_of_pressure(self): + """Evaluates the center of pressure of the aerodynamic surface in local + coordinates. + + Returns + ------- + None + """ + + @abstractmethod + def evaluate_lift_coefficient(self): + """Evaluates the lift coefficient curve of the aerodynamic surface. + + Returns + ------- + None + """ + + @abstractmethod + def evaluate_geometrical_parameters(self): + """Evaluates the geometrical parameters of the aerodynamic surface. + + Returns + ------- + None + """ + + @abstractmethod + def info(self): + """Prints and plots summarized information of the aerodynamic surface. + + Returns + ------- + None + """ + + @abstractmethod + def all_info(self): + """Prints and plots all the available information of the aero surface. + + Returns + ------- + None + """ diff --git a/rocketpy/rocket/aero_surface/air_brakes.py b/rocketpy/rocket/aero_surface/air_brakes.py new file mode 100644 index 000000000..ee4830808 --- /dev/null +++ b/rocketpy/rocket/aero_surface/air_brakes.py @@ -0,0 +1,207 @@ +import warnings + +import numpy as np + +from rocketpy.mathutils.function import Function +from rocketpy.plots.aero_surface_plots import _AirBrakesPlots +from rocketpy.prints.aero_surface_prints import _AirBrakesPrints + +from .aero_surface import AeroSurface + + +class AirBrakes(AeroSurface): + """AirBrakes class. Inherits from AeroSurface. + + Attributes + ---------- + AirBrakes.drag_coefficient : Function + Drag coefficient as a function of deployment level and Mach number. + AirBrakes.drag_coefficient_curve : int, float, callable, array, string, Function + Curve that defines the drag coefficient as a function of deployment level + and Mach number. Used as the source of `AirBrakes.drag_coefficient`. + AirBrakes.deployment_level : float + Current deployment level, ranging from 0 to 1. Deployment level is the + fraction of the total airbrake area that is deployed. + AirBrakes.reference_area : int, float + Reference area used to calculate the drag force of the air brakes + from the drag coefficient curve. Units of m^2. + AirBrakes.clamp : bool, optional + If True, the simulation will clamp the deployment level to 0 or 1 if + the deployment level is out of bounds. If False, the simulation will + not clamp the deployment level and will instead raise a warning if + the deployment level is out of bounds. Default is True. + AirBrakes.name : str + Name of the air brakes. + """ + + def __init__( + self, + drag_coefficient_curve, + reference_area, + clamp=True, + override_rocket_drag=False, + deployment_level=0, + name="AirBrakes", + ): + """Initializes the AirBrakes class. + + Parameters + ---------- + drag_coefficient_curve : int, float, callable, array, string, Function + This parameter represents the drag coefficient associated with the + air brakes and/or the entire rocket, depending on the value of + ``override_rocket_drag``. + + - If a constant, it should be an integer or a float representing a + fixed drag coefficient value. + - If a function, it must take two parameters: deployment level and + Mach number, and return the drag coefficient. This function allows + for dynamic computation based on deployment and Mach number. + - If an array, it should be a 2D array with three columns: the first + column for deployment level, the second for Mach number, and the + third for the corresponding drag coefficient. + - If a string, it should be the path to a .csv or .txt file. The + file must contain three columns: the first for deployment level, + the second for Mach number, and the third for the drag + coefficient. + - If a Function, it must take two parameters: deployment level and + Mach number, and return the drag coefficient. + + .. note:: For ``override_rocket_drag = False``, at + deployment level 0, the drag coefficient is assumed to be 0, + independent of the input drag coefficient curve. This means that + the simulation always considers that at a deployment level of 0, + the air brakes are completely retracted and do not contribute to + the drag of the rocket. + + reference_area : int, float + Reference area used to calculate the drag force of the air brakes + from the drag coefficient curve. Units of m^2. + clamp : bool, optional + If True, the simulation will clamp the deployment level to 0 or 1 if + the deployment level is out of bounds. If False, the simulation will + not clamp the deployment level and will instead raise a warning if + the deployment level is out of bounds. Default is True. + override_rocket_drag : bool, optional + If False, the air brakes drag coefficient will be added to the + rocket's power off drag coefficient curve. If True, during the + simulation, the rocket's power off drag will be ignored and the air + brakes drag coefficient will be used for the entire rocket instead. + Default is False. + deployment_level : float, optional + Initial deployment level, ranging from 0 to 1. Deployment level is + the fraction of the total airbrake area that is Deployment. Default + is 0. + name : str, optional + Name of the air brakes. Default is "AirBrakes". + + Returns + ------- + None + """ + super().__init__(name, reference_area, None) + self.drag_coefficient_curve = drag_coefficient_curve + self.drag_coefficient = Function( + drag_coefficient_curve, + inputs=["Deployment Level", "Mach"], + outputs="Drag Coefficient", + ) + self.clamp = clamp + self.override_rocket_drag = override_rocket_drag + self.initial_deployment_level = deployment_level + self.deployment_level = deployment_level + self.prints = _AirBrakesPrints(self) + self.plots = _AirBrakesPlots(self) + + @property + def deployment_level(self): + """Returns the deployment level of the air brakes.""" + return self._deployment_level + + @deployment_level.setter + def deployment_level(self, value): + # Check if deployment level is within bounds and warn user if not + if value < 0 or value > 1: + # Clamp deployment level if clamp is True + if self.clamp: + # Make sure deployment level is between 0 and 1 + value = np.clip(value, 0, 1) + else: + # Raise warning if clamp is False + warnings.warn( + f"Deployment level of {self.name} is smaller than 0 or " + + "larger than 1. Extrapolation for the drag coefficient " + + "curve will be used." + ) + self._deployment_level = value + + def _reset(self): + """Resets the air brakes to their initial state. This is ran at the + beginning of each simulation to ensure the air brakes are in the correct + state.""" + self.deployment_level = self.initial_deployment_level + + def evaluate_center_of_pressure(self): + """Evaluates the center of pressure of the aerodynamic surface in local + coordinates. + + For air brakes, all components of the center of pressure position are + 0. + + Returns + ------- + None + """ + self.cpx = 0 + self.cpy = 0 + self.cpz = 0 + self.cp = (self.cpx, self.cpy, self.cpz) + + def evaluate_lift_coefficient(self): + """Evaluates the lift coefficient curve of the aerodynamic surface. + + For air brakes, the current model assumes no lift is generated. + Therefore, the lift coefficient (C_L) and its derivative relative to the + angle of attack (C_L_alpha), is 0. + + Returns + ------- + None + """ + self.clalpha = Function( + lambda mach: 0, + "Mach", + f"Lift coefficient derivative for {self.name}", + ) + self.cl = Function( + lambda alpha, mach: 0, + ["Alpha (rad)", "Mach"], + "Lift Coefficient", + ) + + def evaluate_geometrical_parameters(self): + """Evaluates the geometrical parameters of the aerodynamic surface. + + Returns + ------- + None + """ + + def info(self): + """Prints and plots summarized information of the aerodynamic surface. + + Returns + ------- + None + """ + self.prints.geometry() + + def all_info(self): + """Prints and plots all information of the aerodynamic surface. + + Returns + ------- + None + """ + self.info() + self.plots.drag_coefficient_curve() diff --git a/rocketpy/rocket/aero_surface/fins/__init__.py b/rocketpy/rocket/aero_surface/fins/__init__.py new file mode 100644 index 000000000..f1efc603a --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/__init__.py @@ -0,0 +1,3 @@ +from rocketpy.rocket.aero_surface.fins.elliptical_fins import EllipticalFins +from rocketpy.rocket.aero_surface.fins.fins import Fins +from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins diff --git a/rocketpy/rocket/aero_surface/fins/elliptical_fins.py b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py new file mode 100644 index 000000000..5622900d6 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py @@ -0,0 +1,316 @@ +import numpy as np + +from rocketpy.plots.aero_surface_plots import _EllipticalFinsPlots +from rocketpy.prints.aero_surface_prints import _EllipticalFinsPrints + +from .fins import Fins + + +class EllipticalFins(Fins): + """Class that defines and holds information for an elliptical fin set. + + This class inherits from the Fins class. + + Note + ---- + Local coordinate system: Z axis along the longitudinal axis of symmetry, + positive downwards (top -> bottom). Origin located at the top of the root + chord. + + See Also + -------- + Fins + + Attributes + ---------- + EllipticalFins.n : int + Number of fins in fin set. + EllipticalFins.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + EllipticalFins.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees) + EllipticalFins.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + EllipticalFins.changing_attribute_dict : dict + Dictionary that stores the name and the values of the attributes that + may be changed during a simulation. Useful for control systems. + EllipticalFins.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + EllipticalFins.root_chord : float + Fin root chord in meters. + EllipticalFins.span : float + Fin span in meters. + EllipticalFins.name : string + Name of fin set. + EllipticalFins.sweep_length : float + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + EllipticalFins.sweep_angle : float + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. + EllipticalFins.d : float + Reference diameter of the rocket, in meters. + EllipticalFins.ref_area : float + Reference area of the rocket. + EllipticalFins.Af : float + Area of the longitudinal section of each fin in the set. + EllipticalFins.AR : float + Aspect ratio of each fin in the set. + EllipticalFins.gamma_c : float + Fin mid-chord sweep angle. + EllipticalFins.Yma : float + Span wise position of the mean aerodynamic chord. + EllipticalFins.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + EllipticalFins.tau : float + Geometrical relation used to simplify lift and roll calculations. + EllipticalFins.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + EllipticalFins.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + EllipticalFins.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + EllipticalFins.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + EllipticalFins.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + EllipticalFins.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + EllipticalFins.clalpha : float + Lift coefficient slope. Has units of 1/rad. + """ + + def __init__( + self, + n, + root_chord, + span, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Fins", + ): + """Initialize EllipticalFins class. + + Parameters + ---------- + n : int + Number of fins, from 2 to infinity. + root_chord : int, float + Fin root chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference radius to calculate lift coefficient, in meters. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + sweep_length : int, float, optional + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading + edge measured parallel to the rocket centerline. If not given, the + sweep length is assumed to be equal the root chord minus the tip + chord, in which case the fin is a right trapezoid with its base + perpendicular to the rocket's axis. Cannot be used in conjunction + with sweep_angle. + sweep_angle : int, float, optional + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. If not given, the sweep angle is automatically + calculated, in which case the fin is assumed to be a right trapezoid + with its base perpendicular to the rocket's axis. + Cannot be used in conjunction with sweep_length. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of fin set. + + Returns + ------- + None + """ + + super().__init__( + n, + root_chord, + span, + rocket_radius, + cant_angle, + airfoil, + name, + ) + + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + self.prints = _EllipticalFinsPrints(self) + self.plots = _EllipticalFinsPlots(self) + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the fin set in local + coordinates. The center of pressure position is saved and stored as a + tuple. + + Returns + ------- + None + """ + # Center of pressure position in local coordinates + cpz = 0.288 * self.root_chord + self.cpx = 0 + self.cpy = 0 + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements + """Calculates and saves fin set's geometrical parameters such as the + fins' area, aspect ratio and parameters for roll movement. + + Returns + ------- + None + """ + + # Compute auxiliary geometrical parameters + # pylint: disable=invalid-name + Af = (np.pi * self.root_chord / 2 * self.span) / 2 # Fin area + gamma_c = 0 # Zero for elliptical fins + AR = 2 * self.span**2 / Af # Fin aspect ratio + Yma = ( + self.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) + ) # Span wise coord of mean aero chord + roll_geometrical_constant = ( + self.root_chord + * self.span + * ( + 3 * np.pi * self.span**2 + + 32 * self.rocket_radius * self.span + + 12 * np.pi * self.rocket_radius**2 + ) + / 48 + ) + + # Fin–body interference correction parameters + tau = (self.span + self.rocket_radius) / self.rocket_radius + lift_interference_factor = 1 + 1 / tau + if self.span > self.rocket_radius: + roll_damping_interference_factor = 1 + ( + (self.rocket_radius**2) + * ( + 2 + * (self.rocket_radius**2) + * np.sqrt(self.span**2 - self.rocket_radius**2) + * np.log( + ( + 2 + * self.span + * np.sqrt(self.span**2 - self.rocket_radius**2) + + 2 * self.span**2 + ) + / self.rocket_radius + ) + - 2 + * (self.rocket_radius**2) + * np.sqrt(self.span**2 - self.rocket_radius**2) + * np.log(2 * self.span) + + 2 * self.span**3 + - np.pi * self.rocket_radius * self.span**2 + - 2 * (self.rocket_radius**2) * self.span + + np.pi * self.rocket_radius**3 + ) + ) / ( + 2 + * (self.span**2) + * (self.span / 3 + np.pi * self.rocket_radius / 4) + * (self.span**2 - self.rocket_radius**2) + ) + elif self.span < self.rocket_radius: + roll_damping_interference_factor = 1 - ( + self.rocket_radius**2 + * ( + 2 * self.span**3 + - np.pi * self.span**2 * self.rocket_radius + - 2 * self.span * self.rocket_radius**2 + + np.pi * self.rocket_radius**3 + + 2 + * self.rocket_radius**2 + * np.sqrt(-self.span**2 + self.rocket_radius**2) + * np.arctan( + (self.span) / (np.sqrt(-self.span**2 + self.rocket_radius**2)) + ) + - np.pi + * self.rocket_radius**2 + * np.sqrt(-self.span**2 + self.rocket_radius**2) + ) + ) / ( + 2 + * self.span + * (-self.span**2 + self.rocket_radius**2) + * (self.span**2 / 3 + np.pi * self.span * self.rocket_radius / 4) + ) + else: + roll_damping_interference_factor = (28 - 3 * np.pi) / (4 + 3 * np.pi) + + roll_forcing_interference_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2) + / (tau**2 * (tau - 1) ** 2) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + # Store values + # pylint: disable=invalid-name + self.Af = Af # Fin area + self.AR = AR # Fin aspect ratio + self.gamma_c = gamma_c # Mid chord angle + self.Yma = Yma # Span wise coord of mean aero chord + self.roll_geometrical_constant = roll_geometrical_constant + self.tau = tau + self.lift_interference_factor = lift_interference_factor + self.roll_damping_interference_factor = roll_damping_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + + self.evaluate_shape() + + def evaluate_shape(self): + angles = np.arange(0, 180, 5) + x_array = self.root_chord / 2 + self.root_chord / 2 * np.cos(np.radians(angles)) + y_array = self.span * np.sin(np.radians(angles)) + self.shape_vec = [x_array, y_array] + + def info(self): + self.prints.geometry() + self.prints.lift() + + def all_info(self): + self.prints.all() + self.plots.all() diff --git a/rocketpy/rocket/aero_surface/fins/fins.py b/rocketpy/rocket/aero_surface/fins/fins.py new file mode 100644 index 000000000..2de0176b0 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/fins.py @@ -0,0 +1,375 @@ +import numpy as np + +from rocketpy.mathutils.function import Function + +from ..aero_surface import AeroSurface + + +class Fins(AeroSurface): + """Abstract class that holds common methods for the fin classes. + Cannot be instantiated. + + Note + ---- + Local coordinate system: Z axis along the longitudinal axis of symmetry, + positive downwards (top -> bottom). Origin located at the top of the root + chord. + + Attributes + ---------- + Fins.n : int + Number of fins in fin set. + Fins.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, + in meters. + Fins.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + Fins.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + Fins.changing_attribute_dict : dict + Dictionary that stores the name and the values of the attributes that + may be changed during a simulation. Useful for control systems. + Fins.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + Fins.root_chord : float + Fin root chord in meters. + Fins.tip_chord : float + Fin tip chord in meters. + Fins.span : float + Fin span in meters. + Fins.name : string + Name of fin set. + Fins.sweep_length : float + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + Fins.sweep_angle : float + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. + Fins.d : float + Reference diameter of the rocket. Has units of length and is given + in meters. + Fins.ref_area : float + Reference area of the rocket. + Fins.Af : float + Area of the longitudinal section of each fin in the set. + Fins.AR : float + Aspect ratio of each fin in the set. + Fins.gamma_c : float + Fin mid-chord sweep angle. + Fins.Yma : float + Span wise position of the mean aerodynamic chord. + Fins.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + Fins.tau : float + Geometrical relation used to simplify lift and roll calculations. + Fins.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + Fins.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + Fins.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + Fins.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + Fins.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + Fins.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + Fins.clalpha : float + Lift coefficient slope. Has units of 1/rad. + Fins.roll_parameters : list + List containing the roll moment lift coefficient, the roll moment + damping coefficient and the cant angle in radians. + """ + + def __init__( + self, + n, + root_chord, + span, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Fins", + ): + """Initialize Fins class. + + Parameters + ---------- + n : int + Number of fins, from 2 to infinity. + root_chord : int, float + Fin root chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference rocket radius used for lift coefficient normalization. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of fin set. + + Returns + ------- + None + """ + # Compute auxiliary geometrical parameters + d = 2 * rocket_radius + ref_area = np.pi * rocket_radius**2 # Reference area + + super().__init__(name, ref_area, d) + + # Store values + self._n = n + self._rocket_radius = rocket_radius + self._airfoil = airfoil + self._cant_angle = cant_angle + self._root_chord = root_chord + self._span = span + self.name = name + self.d = d + self.ref_area = ref_area # Reference area + + @property + def n(self): + return self._n + + @n.setter + def n(self, value): + self._n = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def root_chord(self): + return self._root_chord + + @root_chord.setter + def root_chord(self, value): + self._root_chord = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def span(self): + return self._span + + @span.setter + def span(self, value): + self._span = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def rocket_radius(self): + return self._rocket_radius + + @rocket_radius.setter + def rocket_radius(self, value): + self._rocket_radius = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def cant_angle(self): + return self._cant_angle + + @cant_angle.setter + def cant_angle(self, value): + self._cant_angle = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def airfoil(self): + return self._airfoil + + @airfoil.setter + def airfoil(self, value): + self._airfoil = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + def evaluate_lift_coefficient(self): + """Calculates and returns the fin set's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves the lift coefficient derivative + for a single fin and the lift coefficient derivative for + a number of n fins corrected for Fin-Body interference. + + Returns + ------- + None + """ + if not self.airfoil: + # Defines clalpha2D as 2*pi for planar fins + clalpha2D_incompressible = 2 * np.pi + else: + # Defines clalpha2D as the derivative of the lift coefficient curve + # for the specific airfoil + self.airfoil_cl = Function( + self.airfoil[0], + interpolation="linear", + ) + + # Differentiating at alpha = 0 to get cl_alpha + clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step( + x=1e-3, dx=1e-3 + ) + + # Convert to radians if needed + if self.airfoil[1] == "degrees": + clalpha2D_incompressible *= 180 / np.pi + + # Correcting for compressible flow (apply Prandtl-Glauert correction) + clalpha2D = Function(lambda mach: clalpha2D_incompressible / self._beta(mach)) + + # Diederich's Planform Correlation Parameter + planform_correlation_parameter = ( + 2 * np.pi * self.AR / (clalpha2D * np.cos(self.gamma_c)) + ) + + # Lift coefficient derivative for a single fin + def lift_source(mach): + return ( + clalpha2D(mach) + * planform_correlation_parameter(mach) + * (self.Af / self.ref_area) + * np.cos(self.gamma_c) + ) / ( + 2 + + planform_correlation_parameter(mach) + * np.sqrt(1 + (2 / planform_correlation_parameter(mach)) ** 2) + ) + + self.clalpha_single_fin = Function( + lift_source, + "Mach", + "Lift coefficient derivative for a single fin", + ) + + # Lift coefficient derivative for n fins corrected with Fin-Body interference + self.clalpha_multiple_fins = ( + self.lift_interference_factor + * self.fin_num_correction(self.n) + * self.clalpha_single_fin + ) # Function of mach number + self.clalpha_multiple_fins.set_inputs("Mach") + self.clalpha_multiple_fins.set_outputs( + f"Lift coefficient derivative for {self.n:.0f} fins" + ) + self.clalpha = self.clalpha_multiple_fins + + # Cl = clalpha * alpha + self.cl = Function( + lambda alpha, mach: alpha * self.clalpha_multiple_fins(mach), + ["Alpha (rad)", "Mach"], + "Lift coefficient", + ) + + return self.cl + + def evaluate_roll_parameters(self): + """Calculates and returns the fin set's roll coefficients. + The roll coefficients are saved in a list. + + Returns + ------- + self.roll_parameters : list + List containing the roll moment lift coefficient, the + roll moment damping coefficient and the cant angle in + radians + """ + + self.cant_angle_rad = np.radians(self.cant_angle) + + clf_delta = ( + self.roll_forcing_interference_factor + * self.n + * (self.Yma + self.rocket_radius) + * self.clalpha_single_fin + / self.d + ) # Function of mach number + clf_delta.set_inputs("Mach") + clf_delta.set_outputs("Roll moment forcing coefficient derivative") + cld_omega = ( + 2 + * self.roll_damping_interference_factor + * self.n + * self.clalpha_single_fin + * np.cos(self.cant_angle_rad) + * self.roll_geometrical_constant + / (self.ref_area * self.d**2) + ) # Function of mach number + cld_omega.set_inputs("Mach") + cld_omega.set_outputs("Roll moment damping coefficient derivative") + self.roll_parameters = [clf_delta, cld_omega, self.cant_angle_rad] + return self.roll_parameters + + @staticmethod + def fin_num_correction(n): + """Calculates a correction factor for the lift coefficient of multiple + fins. + The specifics values are documented at: + Niskanen, S. (2013). “OpenRocket technical documentation”. + In: Development of an Open Source model rocket simulation software. + + Parameters + ---------- + n : int + Number of fins. + + Returns + ------- + Corrector factor : int + Factor that accounts for the number of fins. + """ + corrector_factor = [2.37, 2.74, 2.99, 3.24] + if 5 <= n <= 8: + return corrector_factor[n - 5] + else: + return n / 2 + + def draw(self): + """Draw the fin shape along with some important information, including + the center line, the quarter line and the center of pressure position. + + Returns + ------- + None + """ + self.plots.draw() diff --git a/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py new file mode 100644 index 000000000..15caa38e8 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py @@ -0,0 +1,347 @@ +import numpy as np + +from rocketpy.plots.aero_surface_plots import _TrapezoidalFinsPlots +from rocketpy.prints.aero_surface_prints import _TrapezoidalFinsPrints + +from .fins import Fins + + +class TrapezoidalFins(Fins): + """Class that defines and holds information for a trapezoidal fin set. + + This class inherits from the Fins class. + + Note + ---- + Local coordinate system: Z axis along the longitudinal axis of symmetry, + positive downwards (top -> bottom). Origin located at the top of the root + chord. + + See Also + -------- + Fins + + Attributes + ---------- + TrapezoidalFins.n : int + Number of fins in fin set. + TrapezoidalFins.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + TrapezoidalFins.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + TrapezoidalFins.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + TrapezoidalFins.changing_attribute_dict : dict + Dictionary that stores the name and the values of the attributes that + may be changed during a simulation. Useful for control systems. + TrapezoidalFins.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + TrapezoidalFins.root_chord : float + Fin root chord in meters. + TrapezoidalFins.tip_chord : float + Fin tip chord in meters. + TrapezoidalFins.span : float + Fin span in meters. + TrapezoidalFins.name : string + Name of fin set. + TrapezoidalFins.sweep_length : float + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + TrapezoidalFins.sweep_angle : float + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. + TrapezoidalFins.d : float + Reference diameter of the rocket, in meters. + TrapezoidalFins.ref_area : float + Reference area of the rocket, in m². + TrapezoidalFins.Af : float + Area of the longitudinal section of each fin in the set. + TrapezoidalFins.AR : float + Aspect ratio of each fin in the set + TrapezoidalFins.gamma_c : float + Fin mid-chord sweep angle. + TrapezoidalFins.Yma : float + Span wise position of the mean aerodynamic chord. + TrapezoidalFins.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + TrapezoidalFins.tau : float + Geometrical relation used to simplify lift and roll calculations. + TrapezoidalFins.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + TrapezoidalFins.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + TrapezoidalFins.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + TrapezoidalFins.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + TrapezoidalFins.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + TrapezoidalFins.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + TrapezoidalFins.clalpha : float + Lift coefficient slope. Has units of 1/rad. + """ + + def __init__( + self, + n, + root_chord, + tip_chord, + span, + rocket_radius, + cant_angle=0, + sweep_length=None, + sweep_angle=None, + airfoil=None, + name="Fins", + ): + """Initialize TrapezoidalFins class. + + Parameters + ---------- + n : int + Number of fins, from 2 to infinity. + root_chord : int, float + Fin root chord in meters. + tip_chord : int, float + Fin tip chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference radius to calculate lift coefficient, in meters. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + sweep_length : int, float, optional + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading + edge measured parallel to the rocket centerline. If not given, the + sweep length is assumed to be equal the root chord minus the tip + chord, in which case the fin is a right trapezoid with its base + perpendicular to the rocket's axis. Cannot be used in conjunction + with sweep_angle. + sweep_angle : int, float, optional + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. If not given, the sweep angle is automatically + calculated, in which case the fin is assumed to be a right trapezoid + with its base perpendicular to the rocket's axis. + Cannot be used in conjunction with sweep_length. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of fin set. + + Returns + ------- + None + """ + + super().__init__( + n, + root_chord, + span, + rocket_radius, + cant_angle, + airfoil, + name, + ) + + # Check if sweep angle or sweep length is given + if sweep_length is not None and sweep_angle is not None: + raise ValueError("Cannot use sweep_length and sweep_angle together") + elif sweep_angle is not None: + sweep_length = np.tan(sweep_angle * np.pi / 180) * span + elif sweep_length is None: + sweep_length = root_chord - tip_chord + else: + # Sweep length is given + pass + + self._tip_chord = tip_chord + self._sweep_length = sweep_length + self._sweep_angle = sweep_angle + + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + self.prints = _TrapezoidalFinsPrints(self) + self.plots = _TrapezoidalFinsPlots(self) + + @property + def tip_chord(self): + return self._tip_chord + + @tip_chord.setter + def tip_chord(self, value): + self._tip_chord = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def sweep_angle(self): + return self._sweep_angle + + @sweep_angle.setter + def sweep_angle(self, value): + self._sweep_angle = value + self._sweep_length = np.tan(value * np.pi / 180) * self.span + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def sweep_length(self): + return self._sweep_length + + @sweep_length.setter + def sweep_length(self, value): + self._sweep_length = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the fin set in local + coordinates. The center of pressure position is saved and stored as a + tuple. + + Returns + ------- + None + """ + # Center of pressure position in local coordinates + cpz = (self.sweep_length / 3) * ( + (self.root_chord + 2 * self.tip_chord) / (self.root_chord + self.tip_chord) + ) + (1 / 6) * ( + self.root_chord + + self.tip_chord + - self.root_chord * self.tip_chord / (self.root_chord + self.tip_chord) + ) + self.cpx = 0 + self.cpy = 0 + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements + """Calculates and saves fin set's geometrical parameters such as the + fins' area, aspect ratio and parameters for roll movement. + + Returns + ------- + None + """ + # pylint: disable=invalid-name + Yr = self.root_chord + self.tip_chord + Af = Yr * self.span / 2 # Fin area + AR = 2 * self.span**2 / Af # Fin aspect ratio + gamma_c = np.arctan( + (self.sweep_length + 0.5 * self.tip_chord - 0.5 * self.root_chord) + / (self.span) + ) + Yma = ( + (self.span / 3) * (self.root_chord + 2 * self.tip_chord) / Yr + ) # Span wise coord of mean aero chord + + # Fin–body interference correction parameters + tau = (self.span + self.rocket_radius) / self.rocket_radius + lift_interference_factor = 1 + 1 / tau + lambda_ = self.tip_chord / self.root_chord + + # Parameters for Roll Moment. + # Documented at: https://docs.rocketpy.org/en/latest/technical/ + roll_geometrical_constant = ( + (self.root_chord + 3 * self.tip_chord) * self.span**3 + + 4 + * (self.root_chord + 2 * self.tip_chord) + * self.rocket_radius + * self.span**2 + + 6 * (self.root_chord + self.tip_chord) * self.span * self.rocket_radius**2 + ) / 12 + roll_damping_interference_factor = 1 + ( + ((tau - lambda_) / (tau)) - ((1 - lambda_) / (tau - 1)) * np.log(tau) + ) / ( + ((tau + 1) * (tau - lambda_)) / (2) + - ((1 - lambda_) * (tau**3 - 1)) / (3 * (tau - 1)) + ) + roll_forcing_interference_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2) + / (tau**2 * (tau - 1) ** 2) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + # Store values + self.Yr = Yr + self.Af = Af # Fin area + self.AR = AR # Aspect Ratio + self.gamma_c = gamma_c # Mid chord angle + self.Yma = Yma # Span wise coord of mean aero chord + self.roll_geometrical_constant = roll_geometrical_constant + self.tau = tau + self.lift_interference_factor = lift_interference_factor + self.λ = lambda_ # pylint: disable=non-ascii-name + self.roll_damping_interference_factor = roll_damping_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + + self.evaluate_shape() + + def evaluate_shape(self): + if self.sweep_length: + points = [ + (0, 0), + (self.sweep_length, self.span), + (self.sweep_length + self.tip_chord, self.span), + (self.root_chord, 0), + ] + else: + points = [ + (0, 0), + (self.root_chord - self.tip_chord, self.span), + (self.root_chord, self.span), + (self.root_chord, 0), + ] + + x_array, y_array = zip(*points) + self.shape_vec = [np.array(x_array), np.array(y_array)] + + def info(self): + self.prints.geometry() + self.prints.lift() + + def all_info(self): + self.prints.all() + self.plots.all() diff --git a/rocketpy/rocket/aero_surface/nose_cone.py b/rocketpy/rocket/aero_surface/nose_cone.py new file mode 100644 index 000000000..8886be8c5 --- /dev/null +++ b/rocketpy/rocket/aero_surface/nose_cone.py @@ -0,0 +1,521 @@ +import warnings + +import numpy as np +from scipy.optimize import fsolve + +from rocketpy.mathutils.function import Function +from rocketpy.plots.aero_surface_plots import _NoseConePlots +from rocketpy.prints.aero_surface_prints import _NoseConePrints + +from .aero_surface import AeroSurface + + +class NoseCone(AeroSurface): + """Keeps nose cone information. + + Note + ---- + The local coordinate system has the origin at the tip of the nose cone + and the Z axis along the longitudinal axis of symmetry, positive + downwards (top -> bottom). + + Attributes + ---------- + NoseCone.length : float + Nose cone length. Has units of length and must be given in meters. + NoseCone.kind : string + Nose cone kind. Can be "conical", "ogive", "elliptical", "tangent", + "von karman", "parabolic", "powerseries" or "lvhaack". + NoseCone.bluffness : float + Ratio between the radius of the circle on the tip of the ogive and the + radius of the base of the ogive. Currently only used for the nose cone's + drawing. Must be between 0 and 1. Default is None, which means that the + nose cone will not have a sphere on the tip. If a value is given, the + nose cone's length will be slightly reduced because of the addition of + the sphere. Must be None or 0 if a "powerseries" nose cone kind is + specified. + NoseCone.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, + in meters. + NoseCone.base_radius : float + Nose cone base radius. Has units of length and must be given in meters. + NoseCone.radius_ratio : float + Ratio between the nose cone base radius and the rocket radius. Has no + units. If base radius is not given, the ratio between base radius and + rocket radius is assumed as 1, meaning that the nose cone has the same + radius as the rocket. If base radius is given, the ratio between base + radius and rocket radius is calculated and used for lift calculation. + NoseCone.power : float + Factor that controls the bluntness of the shape for a power series + nose cone. Must be between 0 and 1. It is ignored when other nose + cone types are used. + NoseCone.name : string + Nose cone name. Has no impact in simulation, as it is only used to + display data in a more organized matter. + NoseCone.cp : tuple + Tuple with the x, y and z local coordinates of the nose cone center of + pressure. Has units of length and is given in meters. + NoseCone.cpx : float + Nose cone local center of pressure x coordinate. Has units of length and + is given in meters. + NoseCone.cpy : float + Nose cone local center of pressure y coordinate. Has units of length and + is given in meters. + NoseCone.cpz : float + Nose cone local center of pressure z coordinate. Has units of length and + is given in meters. + NoseCone.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + NoseCone.clalpha : float + Lift coefficient slope. Has units of 1/rad. + NoseCone.plots : plots.aero_surface_plots._NoseConePlots + This contains all the plots methods. Use help(NoseCone.plots) to know + more about it. + NoseCone.prints : prints.aero_surface_prints._NoseConePrints + This contains all the prints methods. Use help(NoseCone.prints) to know + more about it. + """ + + def __init__( # pylint: disable=too-many-statements + self, + length, + kind, + base_radius=None, + bluffness=None, + rocket_radius=None, + power=None, + name="Nose Cone", + ): + """Initializes the nose cone. It is used to define the nose cone + length, kind, center of pressure and lift coefficient curve. + + Parameters + ---------- + length : float + Nose cone length. Has units of length and must be given in meters. + kind : string + Nose cone kind. Can be "conical", "ogive", "elliptical", "tangent", + "von karman", "parabolic", "powerseries" or "lvhaack". If + "powerseries" is used, the "power" argument must be assigned to a + value between 0 and 1. + base_radius : float, optional + Nose cone base radius. Has units of length and must be given in + meters. If not given, the ratio between ``base_radius`` and + ``rocket_radius`` will be assumed as 1. + bluffness : float, optional + Ratio between the radius of the circle on the tip of the ogive and + the radius of the base of the ogive. Currently only used for the + nose cone's drawing. Must be between 0 and 1. Default is None, which + means that the nose cone will not have a sphere on the tip. If a + value is given, the nose cone's length will be reduced to account + for the addition of the sphere at the tip. Must be None or 0 if a + "powerseries" nose cone kind is specified. + rocket_radius : int, float, optional + The reference rocket radius used for lift coefficient normalization. + If not given, the ratio between ``base_radius`` and + ``rocket_radius`` will be assumed as 1. + power : float, optional + Factor that controls the bluntness of the shape for a power series + nose cone. Must be between 0 and 1. It is ignored when other nose + cone types are used. + name : str, optional + Nose cone name. Has no impact in simulation, as it is only used to + display data in a more organized matter. + + Returns + ------- + None + """ + rocket_radius = rocket_radius or base_radius + super().__init__(name, np.pi * rocket_radius**2, 2 * rocket_radius) + + self._rocket_radius = rocket_radius + self._base_radius = base_radius + self._length = length + if bluffness is not None: + if bluffness > 1 or bluffness < 0: + raise ValueError( + f"Bluffness ratio of {bluffness} is out of range. " + "It must be between 0 and 1." + ) + self._bluffness = bluffness + if kind == "powerseries": + # Checks if bluffness is not being used + if (self.bluffness is not None) and (self.bluffness != 0): + raise ValueError( + "Parameter 'bluffness' must be None or 0 when using a nose cone kind 'powerseries'." + ) + + if power is None: + raise ValueError( + "Parameter 'power' cannot be None when using a nose cone kind 'powerseries'." + ) + + if power > 1 or power <= 0: + raise ValueError( + f"Power value of {power} is out of range. It must be between 0 and 1." + ) + self._power = power + self.kind = kind + + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + + self.plots = _NoseConePlots(self) + self.prints = _NoseConePrints(self) + + @property + def rocket_radius(self): + return self._rocket_radius + + @rocket_radius.setter + def rocket_radius(self, value): + self._rocket_radius = value + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_nose_shape() + + @property + def base_radius(self): + return self._base_radius + + @base_radius.setter + def base_radius(self, value): + self._base_radius = value + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_nose_shape() + + @property + def length(self): + return self._length + + @length.setter + def length(self, value): + self._length = value + self.evaluate_center_of_pressure() + self.evaluate_nose_shape() + + @property + def power(self): + return self._power + + @power.setter + def power(self, value): + if value is not None: + if value > 1 or value <= 0: + raise ValueError( + f"Power value of {value} is out of range. It must be between 0 and 1." + ) + self._power = value + self.evaluate_k() + self.evaluate_center_of_pressure() + self.evaluate_nose_shape() + + @property + def kind(self): + return self._kind + + @kind.setter + def kind(self, value): # pylint: disable=too-many-statements + # Analyzes nosecone type + # Sets the k for Cp calculation + # Sets the function which creates the respective curve + self._kind = value + value = (value.replace(" ", "")).lower() + + if value == "conical": + self.k = 2 / 3 + self.y_nosecone = Function(lambda x: x * self.base_radius / self.length) + + elif value == "lvhaack": + self.k = 0.563 + + def theta(x): + return np.arccos(1 - 2 * max(min(x / self.length, 1), 0)) + + self.y_nosecone = Function( + lambda x: self.base_radius + * (theta(x) - np.sin(2 * theta(x)) / 2 + (np.sin(theta(x)) ** 3) / 3) + ** (0.5) + / (np.pi**0.5) + ) + + elif value in ["tangent", "tangentogive", "ogive"]: + rho = (self.base_radius**2 + self.length**2) / (2 * self.base_radius) + volume = np.pi * ( + self.length * rho**2 + - (self.length**3) / 3 + - (rho - self.base_radius) * rho**2 * np.arcsin(self.length / rho) + ) + area = np.pi * self.base_radius**2 + self.k = 1 - volume / (area * self.length) + self.y_nosecone = Function( + lambda x: np.sqrt(rho**2 - (min(x - self.length, 0)) ** 2) + + (self.base_radius - rho) + ) + + elif value == "elliptical": + self.k = 1 / 3 + self.y_nosecone = Function( + lambda x: self.base_radius + * np.sqrt(1 - ((x - self.length) / self.length) ** 2) + ) + + elif value == "vonkarman": + self.k = 0.5 + + def theta(x): + return np.arccos(1 - 2 * max(min(x / self.length, 1), 0)) + + self.y_nosecone = Function( + lambda x: self.base_radius + * (theta(x) - np.sin(2 * theta(x)) / 2) ** (0.5) + / (np.pi**0.5) + ) + elif value == "parabolic": + self.k = 0.5 + self.y_nosecone = Function( + lambda x: self.base_radius + * ((2 * x / self.length - (x / self.length) ** 2) / (2 - 1)) + ) + elif value == "powerseries": + self.k = (2 * self.power) / ((2 * self.power) + 1) + self.y_nosecone = Function( + lambda x: self.base_radius * np.power(x / self.length, self.power) + ) + else: + raise ValueError( + f"Nose Cone kind '{self.kind}' not found, " + + "please use one of the following Nose Cone kinds:" + + '\n\t"conical"' + + '\n\t"ogive"' + + '\n\t"lvhaack"' + + '\n\t"tangent"' + + '\n\t"vonkarman"' + + '\n\t"elliptical"' + + '\n\t"powerseries"' + + '\n\t"parabolic"\n' + ) + + self.evaluate_center_of_pressure() + self.evaluate_geometrical_parameters() + self.evaluate_nose_shape() + + @property + def bluffness(self): + return self._bluffness + + @bluffness.setter + def bluffness(self, value): + # prevents from setting bluffness on "powerseries" nose cones + if self.kind == "powerseries": + # Checks if bluffness is not being used + if (value is not None) and (value != 0): + raise ValueError( + "Parameter 'bluffness' must be None or 0 when using a nose cone kind 'powerseries'." + ) + if value is not None: + if value > 1 or value < 0: + raise ValueError( + f"Bluffness ratio of {value} is out of range. " + "It must be between 0 and 1." + ) + self._bluffness = value + self.evaluate_nose_shape() + + def evaluate_geometrical_parameters(self): + """Calculates and saves nose cone's radius ratio. + + Returns + ------- + None + """ + + # If base radius is not given, the ratio between base radius and + # rocket radius is assumed as 1, meaning that the nose cone has the + # same radius as the rocket + if self.base_radius is None and self.rocket_radius is not None: + self.radius_ratio = 1 + self.base_radius = self.rocket_radius + elif self.base_radius is not None and self.rocket_radius is None: + self.radius_ratio = 1 + self.rocket_radius = self.base_radius + # If base radius is given, the ratio between base radius and rocket + # radius is calculated + elif self.base_radius is not None and self.rocket_radius is not None: + self.radius_ratio = self.base_radius / self.rocket_radius + else: + raise ValueError( + "Either base radius or rocket radius must be given to " + "calculate the nose cone radius ratio." + ) + + self.fineness_ratio = self.length / (2 * self.base_radius) + + def evaluate_nose_shape(self): # pylint: disable=too-many-statements + """Calculates and saves nose cone's shape as lists and re-evaluates the + nose cone's length for a given bluffness ratio. The shape is saved as + two vectors, one for the x coordinates and one for the y coordinates. + + Returns + ------- + None + """ + number_of_points = 127 + density_modifier = 3 # increase density of points to improve accuracy + + def find_x_intercept(x): + # find the tangential intersection point between the circle and nosec curve + return x + self.y_nosecone(x) * self.y_nosecone.differentiate_complex_step( + x + ) + + # Calculate a function to find the radius of the nosecone curve + def find_radius(x): + return (self.y_nosecone(x) ** 2 + (x - find_x_intercept(x)) ** 2) ** 0.5 + + # Check bluffness circle and choose whether to use it or not + if self.bluffness is None or self.bluffness == 0: + # Set up parameters to continue without bluffness + r_circle, circle_center, x_init = 0, 0, 0 + else: + # Calculate circle radius + r_circle = self.bluffness * self.base_radius + if self.kind == "elliptical": + + def test_circle(x): + # set up a circle at the starting position to test bluffness + return np.sqrt(r_circle**2 - (x - r_circle) ** 2) + + # Check if bluffness circle is too small + if test_circle(1e-03) <= self.y_nosecone(1e-03): + # Raise a warning + warnings.warn( + "WARNING: The chosen bluffness ratio is too small for " + "the selected nose cone category, thereby the effective " + "bluffness will be 0." + ) + # Set up parameters to continue without bluffness + r_circle, circle_center, x_init = 0, 0, 0 + else: + # Find the intersection point between circle and nosecone curve + x_init = fsolve(lambda x: find_radius(x[0]) - r_circle, r_circle)[0] + circle_center = find_x_intercept(x_init) + else: + # Find the intersection point between circle and nosecone curve + x_init = fsolve(lambda x: find_radius(x[0]) - r_circle, r_circle)[0] + circle_center = find_x_intercept(x_init) + + # Calculate a function to create the circle at the correct position + def create_circle(x): + return abs(r_circle**2 - (x - circle_center) ** 2) ** 0.5 + + # Define a function for the final shape of the curve with a circle at the tip + def final_shape(x): + return self.y_nosecone(x) if x >= x_init else create_circle(x) + + # Vectorize the final_shape function + final_shape_vec = np.vectorize(final_shape) + + # Create the vectors X and Y with the points of the curve + nosecone_x = (self.length - (circle_center - r_circle)) * ( + np.linspace(0, 1, number_of_points) ** density_modifier + ) + nosecone_y = final_shape_vec(nosecone_x + (circle_center - r_circle)) + + # Evaluate final geometry parameters + self.shape_vec = [nosecone_x, nosecone_y] + if abs(nosecone_x[-1] - self.length) >= 0.001: # 1 millimeter + self._length = nosecone_x[-1] + print( + "Due to the chosen bluffness ratio, the nose " + f"cone length was reduced to {self.length} m." + ) + self.fineness_ratio = self.length / (2 * self.base_radius) + + def evaluate_lift_coefficient(self): + """Calculates and returns nose cone's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves its lift coefficient derivative. + + Returns + ------- + None + """ + # Calculate clalpha + # clalpha is currently a constant, meaning it is independent of Mach + # number. This is only valid for subsonic speeds. + # It must be set as a Function because it will be called and treated + # as a function of mach in the simulation. + self.clalpha = Function( + lambda mach: 2 / self._beta(mach) * self.radius_ratio**2, + "Mach", + f"Lift coefficient derivative for {self.name}", + ) + self.cl = Function( + lambda alpha, mach: self.clalpha(mach) * alpha, + ["Alpha (rad)", "Mach"], + "Cl", + ) + + def evaluate_k(self): + """Updates the self.k attribute used to compute the center of + pressure when using "powerseries" nose cones. + + Returns + ------- + None + """ + if self.kind == "powerseries": + self.k = (2 * self.power) / ((2 * self.power) + 1) + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the nose cone in + local coordinates. The center of pressure position is saved and stored + as a tuple. Local coordinate origin is found at the tip of the nose + cone. + + Returns + ------- + self.cp : tuple + Tuple containing cpx, cpy, cpz. + """ + + self.cpz = self.k * self.length + self.cpy = 0 + self.cpx = 0 + self.cp = (self.cpx, self.cpy, self.cpz) + return self.cp + + def draw(self): + """Draw the fin shape along with some important information, including + the center line, the quarter line and the center of pressure position. + + Returns + ------- + None + """ + self.plots.draw() + + def info(self): + """Prints and plots summarized information of the nose cone. + + Return + ------ + None + """ + self.prints.geometry() + self.prints.lift() + + def all_info(self): + """Prints and plots all the available information of the nose cone. + + Returns + ------- + None + """ + self.prints.all() + self.plots.all() diff --git a/rocketpy/rocket/aero_surface/rail_buttons.py b/rocketpy/rocket/aero_surface/rail_buttons.py new file mode 100644 index 000000000..6ffafbb97 --- /dev/null +++ b/rocketpy/rocket/aero_surface/rail_buttons.py @@ -0,0 +1,106 @@ +from rocketpy.mathutils.function import Function +from rocketpy.prints.aero_surface_prints import _RailButtonsPrints + +from .aero_surface import AeroSurface + + +class RailButtons(AeroSurface): + """Class that defines a rail button pair or group. + + Attributes + ---------- + RailButtons.buttons_distance : int, float + Distance between the two rail buttons closest to the nozzle. + RailButtons.angular_position : int, float + Angular position of the rail buttons in degrees measured + as the rotation around the symmetry axis of the rocket + relative to one of the other principal axis. + """ + + def __init__(self, buttons_distance, angular_position=45, name="Rail Buttons"): + """Initializes RailButtons Class. + + Parameters + ---------- + buttons_distance : int, float + Distance between the first and the last rail button in meters. + angular_position : int, float, optional + Angular position of the rail buttons in degrees measured + as the rotation around the symmetry axis of the rocket + relative to one of the other principal axis. + name : string, optional + Name of the rail buttons. Default is "Rail Buttons". + + Returns + ------- + None + + """ + super().__init__(name, None, None) + self.buttons_distance = buttons_distance + self.angular_position = angular_position + self.name = name + + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + + self.prints = _RailButtonsPrints(self) + + def evaluate_center_of_pressure(self): + """Evaluates the center of pressure of the rail buttons. Rail buttons + do not contribute to the center of pressure of the rocket. + + Returns + ------- + None + """ + self.cpx = 0 + self.cpy = 0 + self.cpz = 0 + self.cp = (self.cpx, self.cpy, self.cpz) + + def evaluate_lift_coefficient(self): + """Evaluates the lift coefficient curve of the rail buttons. Rail + buttons do not contribute to the lift coefficient of the rocket. + + Returns + ------- + None + """ + self.clalpha = Function( + lambda mach: 0, + "Mach", + f"Lift coefficient derivative for {self.name}", + ) + self.cl = Function( + lambda alpha, mach: 0, + ["Alpha (rad)", "Mach"], + "Cl", + ) + + def evaluate_geometrical_parameters(self): + """Evaluates the geometrical parameters of the rail buttons. Rail + buttons do not contribute to the geometrical parameters of the rocket. + + Returns + ------- + None + """ + + def info(self): + """Prints out all the information about the Rail Buttons. + + Returns + ------- + None + """ + self.prints.geometry() + + def all_info(self): + """Returns all info of the Rail Buttons. + + Returns + ------- + None + """ + self.prints.all() diff --git a/rocketpy/rocket/aero_surface/tail.py b/rocketpy/rocket/aero_surface/tail.py new file mode 100644 index 000000000..a5c5cf939 --- /dev/null +++ b/rocketpy/rocket/aero_surface/tail.py @@ -0,0 +1,211 @@ +import numpy as np + +from rocketpy.mathutils.function import Function +from rocketpy.plots.aero_surface_plots import _TailPlots +from rocketpy.prints.aero_surface_prints import _TailPrints + +from .aero_surface import AeroSurface + + +class Tail(AeroSurface): + """Class that defines a tail. Currently only accepts conical tails. + + Note + ---- + Local coordinate system: Z axis along the longitudinal axis of symmetry, + positive downwards (top -> bottom). Origin located at top of the tail + (generally the portion closest to the rocket's nose). + + Attributes + ---------- + Tail.top_radius : int, float + Radius of the top of the tail. The top radius is defined as the radius + of the transversal section that is closest to the rocket's nose. + Tail.bottom_radius : int, float + Radius of the bottom of the tail. + Tail.length : int, float + Length of the tail. The length is defined as the distance between the + top and bottom of the tail. The length is measured along the rocket's + longitudinal axis. Has the unit of meters. + Tail.rocket_radius: int, float + The reference rocket radius used for lift coefficient normalization in + meters. + Tail.name : str + Name of the tail. Default is 'Tail'. + Tail.cpx : int, float + x local coordinate of the center of pressure of the tail. + Tail.cpy : int, float + y local coordinate of the center of pressure of the tail. + Tail.cpz : int, float + z local coordinate of the center of pressure of the tail. + Tail.cp : tuple + Tuple containing the coordinates of the center of pressure of the tail. + Tail.cl : Function + Function that returns the lift coefficient of the tail. The function + is defined as a function of the angle of attack and the mach number. + Tail.clalpha : float + Lift coefficient slope. Has the unit of 1/rad. + Tail.slant_length : float + Slant length of the tail. The slant length is defined as the distance + between the top and bottom of the tail. The slant length is measured + along the tail's slant axis. Has the unit of meters. + Tail.surface_area : float + Surface area of the tail. Has the unit of meters squared. + """ + + def __init__(self, top_radius, bottom_radius, length, rocket_radius, name="Tail"): + """Initializes the tail object by computing and storing the most + important values. + + Parameters + ---------- + top_radius : int, float + Radius of the top of the tail. The top radius is defined as the + radius of the transversal section that is closest to the rocket's + nose. + bottom_radius : int, float + Radius of the bottom of the tail. + length : int, float + Length of the tail. + rocket_radius : int, float + The reference rocket radius used for lift coefficient normalization. + name : str + Name of the tail. Default is 'Tail'. + + Returns + ------- + None + """ + super().__init__(name, np.pi * rocket_radius**2, 2 * rocket_radius) + + self._top_radius = top_radius + self._bottom_radius = bottom_radius + self._length = length + self._rocket_radius = rocket_radius + + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + + self.plots = _TailPlots(self) + self.prints = _TailPrints(self) + + @property + def top_radius(self): + return self._top_radius + + @top_radius.setter + def top_radius(self, value): + self._top_radius = value + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + + @property + def bottom_radius(self): + return self._bottom_radius + + @bottom_radius.setter + def bottom_radius(self, value): + self._bottom_radius = value + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + + @property + def length(self): + return self._length + + @length.setter + def length(self, value): + self._length = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + + @property + def rocket_radius(self): + return self._rocket_radius + + @rocket_radius.setter + def rocket_radius(self, value): + self._rocket_radius = value + self.evaluate_lift_coefficient() + + def evaluate_geometrical_parameters(self): + """Calculates and saves tail's slant length and surface area. + + Returns + ------- + None + """ + self.slant_length = np.sqrt( + (self.length) ** 2 + (self.top_radius - self.bottom_radius) ** 2 + ) + self.surface_area = ( + np.pi * self.slant_length * (self.top_radius + self.bottom_radius) + ) + self.evaluate_shape() + + def evaluate_shape(self): + # Assuming the tail is a cone, calculate the shape vector + self.shape_vec = [ + np.array([0, self.length]), + np.array([self.top_radius, self.bottom_radius]), + ] + + def evaluate_lift_coefficient(self): + """Calculates and returns tail's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves its lift coefficient derivative. + + Returns + ------- + None + """ + # Calculate clalpha + # clalpha is currently a constant, meaning it is independent of Mach + # number. This is only valid for subsonic speeds. + # It must be set as a Function because it will be called and treated + # as a function of mach in the simulation. + self.clalpha = Function( + lambda mach: 2 + / self._beta(mach) + * ( + (self.bottom_radius / self.rocket_radius) ** 2 + - (self.top_radius / self.rocket_radius) ** 2 + ), + "Mach", + f"Lift coefficient derivative for {self.name}", + ) + self.cl = Function( + lambda alpha, mach: self.clalpha(mach) * alpha, + ["Alpha (rad)", "Mach"], + "Cl", + ) + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the tail in local + coordinates. The center of pressure position is saved and stored as a + tuple. + + Returns + ------- + None + """ + # Calculate cp position in local coordinates + r = self.top_radius / self.bottom_radius + cpz = (self.length / 3) * (1 + (1 - r) / (1 - r**2)) + + # Store values as class attributes + self.cpx = 0 + self.cpy = 0 + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + def info(self): + self.prints.geometry() + self.prints.lift() + + def all_info(self): + self.prints.all() + self.plots.all() diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index f5bfc8666..bcaa31cc3 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1440,8 +1440,8 @@ def u_dot( comp_cp = ( position - self.rocket.center_of_dry_mass_position ) * self.rocket._csys - aero_surface.cpz - surface_radius = aero_surface.rocket_radius - reference_area = np.pi * surface_radius**2 + reference_area = aero_surface.reference_area + reference_length = aero_surface.reference_length # Component absolute velocity in body frame comp_vx_b = vx_b + comp_cp * omega2 comp_vy_b = vy_b - comp_cp * omega1 @@ -1491,15 +1491,14 @@ def u_dot( M3_forcing = ( (1 / 2 * rho * free_stream_speed**2) * reference_area - * 2 - * surface_radius + * reference_length * clf_delta.get_value_opt(free_stream_mach) * cant_angle_rad ) M3_damping = ( (1 / 2 * rho * free_stream_speed) * reference_area - * (2 * surface_radius) ** 2 + * (reference_length) ** 2 * cld_omega.get_value_opt(free_stream_mach) * omega3 / 2 @@ -1723,8 +1722,8 @@ def u_dot_generalized( position - self.rocket.center_of_dry_mass_position ) * self.rocket._csys - aero_surface.cpz comp_cp = Vector([0, 0, comp_cpz]) - surface_radius = aero_surface.rocket_radius - reference_area = np.pi * surface_radius**2 + reference_area = aero_surface.reference_area + reference_length = aero_surface.reference_length # Component absolute velocity in body frame comp_vb = velocity_in_body_frame + (w ^ comp_cp) # Wind velocity at component altitude @@ -1768,15 +1767,14 @@ def u_dot_generalized( M3_forcing = ( (1 / 2 * rho * comp_stream_speed**2) * reference_area - * 2 - * surface_radius + * reference_length * clf_delta.get_value_opt(comp_stream_mach) * cant_angle_rad ) M3_damping = ( (1 / 2 * rho * comp_stream_speed) * reference_area - * (2 * surface_radius) ** 2 + * (reference_length) ** 2 * cld_omega.get_value_opt(comp_stream_mach) * omega3 / 2