diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py new file mode 100644 index 000000000..8f64bf1b1 --- /dev/null +++ b/rocketpy/rocket/aero_surface.py @@ -0,0 +1,2152 @@ +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, rocket_radius, 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". + rocket_radius: int, float, optional + The reference rocket radius used for lift coefficient normalization. + + Returns + ------- + None + + """ + super().__init__(name) + self.buttons_distance = buttons_distance + self.angular_position = angular_position + self.name = name + self.rocket_radius = rocket_radius + 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/rocket.py b/rocketpy/rocket/rocket.py index d6dff9113..d64c3badd 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1509,8 +1509,9 @@ def set_rail_buttons( RailButtons object created """ buttons_distance = abs(upper_button_position - lower_button_position) - rail_buttons = RailButtons( - buttons_distance=buttons_distance, angular_position=angular_position + rail_buttons = RailButtons(rocket_radius=self.radius, + buttons_distance=buttons_distance, + angular_position=angular_position ) self.rail_buttons.add(rail_buttons, lower_button_position) return rail_buttons