diff --git a/docs/notebooks/coeff_testing.ipynb b/docs/notebooks/coeff_testing.ipynb index 80ea83b55..5bd822164 100644 --- a/docs/notebooks/coeff_testing.ipynb +++ b/docs/notebooks/coeff_testing.ipynb @@ -458,25 +458,31 @@ "gennose = GenericSurface(\n", " reference_area=np.pi * calisto.radius**2,\n", " reference_length=2 * calisto.radius,\n", - " cL=\"nose_cL.csv\",\n", - " cQ=\"nose_cQ.csv\",\n", + " coefficients={\n", + " \"cL\": \"nose_cL.csv\",\n", + " \"cQ\": \"nose_cQ.csv\",\n", + " },\n", " center_of_pressure=(0, 0, 0),\n", " name=\"nose\",\n", ")\n", "genfin = GenericSurface(\n", " reference_area=np.pi * calisto.radius**2,\n", " reference_length=2 * calisto.radius,\n", - " cL=\"fins_cL.csv\",\n", - " cQ=\"fins_cQ.csv\",\n", - " cl=\"fins_roll.csv\",\n", + " coefficients={\n", + " \"cL\": \"fins_cL.csv\",\n", + " \"cQ\": \"fins_cQ.csv\",\n", + " \"cl\": \"fins_roll.csv\",\n", + " },\n", " center_of_pressure=(0, 0, 0),\n", " name=\"fins\",\n", ")\n", "gentail = GenericSurface(\n", " reference_area=np.pi * calisto.radius**2,\n", " reference_length=2 * calisto.radius,\n", - " cL=\"tail_cL.csv\",\n", - " cQ=\"tail_cQ.csv\",\n", + " coefficients={\n", + " \"cL\": \"tail_cL.csv\",\n", + " \"cQ\": \"tail_cQ.csv\",\n", + " },\n", " center_of_pressure=(0, 0, 0),\n", " name=\"tail\",\n", ")\n", @@ -2451,7 +2457,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.0" + "version": "3.11.2" }, "vscode": { "interpreter": { diff --git a/rocketpy/rocket/aero_surface/aero_surface.py b/rocketpy/rocket/aero_surface/aero_surface.py index 5b98f6260..15ca14f1d 100644 --- a/rocketpy/rocket/aero_surface/aero_surface.py +++ b/rocketpy/rocket/aero_surface/aero_surface.py @@ -99,8 +99,7 @@ def compute_forces_and_moments( rho, cp, *args, - **kwargs, - ): + ): # pylint: disable=unused-argument """Computes the forces and moments acting on the aerodynamic surface. Used in each time step of the simulation. This method is valid for the barrowman aerodynamic models. diff --git a/rocketpy/rocket/aero_surface/fins/fins.py b/rocketpy/rocket/aero_surface/fins/fins.py index b4aa2f790..e6cfb84ea 100644 --- a/rocketpy/rocket/aero_surface/fins/fins.py +++ b/rocketpy/rocket/aero_surface/fins/fins.py @@ -373,21 +373,33 @@ def compute_forces_and_moments( stream_mach, rho, cp, - _, - omega1, - omega2, - omega3, + omega, *args, - **kwargs, - ): + ): # pylint: disable=arguments-differ """Computes the forces and moments acting on the aerodynamic surface. Parameters ---------- - stream_speed : int, float - Speed of the flow stream in the body frame. + stream_velocity : tuple of float + The velocity of the airflow relative to the surface. + stream_speed : float + The magnitude of the airflow speed. + stream_mach : float + The Mach number of the airflow. + rho : float + Air density. + cp : Vector + Center of pressure coordinates in the body frame. + omega: tuple[float, float, float] + Tuple containing angular velocities around the x, y, z axes. + Returns + ------- + tuple of float + The aerodynamic forces (lift, side_force, drag) and moments + (pitch, yaw, roll) in the body frame. """ + R1, R2, R3, M1, M2, _ = super().compute_forces_and_moments( stream_velocity, stream_speed, @@ -408,7 +420,7 @@ def compute_forces_and_moments( * self.reference_area * (self.reference_length) ** 2 * cld_omega.get_value_opt(stream_mach) - * omega3 + * omega[2] / 2 ) M3 = M3_forcing - M3_damping diff --git a/rocketpy/rocket/aero_surface/generic_surface.py b/rocketpy/rocket/aero_surface/generic_surface.py index ff86bd942..b7c803fa3 100644 --- a/rocketpy/rocket/aero_surface/generic_surface.py +++ b/rocketpy/rocket/aero_surface/generic_surface.py @@ -1,3 +1,4 @@ +import copy import csv import math @@ -17,12 +18,7 @@ def __init__( self, reference_area, reference_length, - cL=0, - cQ=0, - cD=0, - cm=0, - cn=0, - cl=0, + coefficients, center_of_pressure=(0, 0, 0), name="Generic Surface", ): @@ -51,24 +47,27 @@ def __init__( reference_length : int, float Reference length of the aerodynamic surface. Has the unit of meters. Commonly defined as the rocket's diameter. - cL : str, callable, optional - Lift coefficient. Can be a path to a CSV file or a callable. - Default is 0. - cQ : str, callable, optional - Side force coefficient. Can be a path to a CSV file or a callable. - Default is 0. - cD : str, callable, optional - Drag coefficient. Can be a path to a CSV file or a callable. - Default is 0. - cm : str, callable, optional - Pitch moment coefficient. Can be a path to a CSV file or a callable. - Default is 0. - cn : str, callable, optional - Yaw moment coefficient. Can be a path to a CSV file or a callable. - Default is 0. - cl : str, callable, optional - Roll moment coefficient. Can be a path to a CSV file or a callable. - Default is 0. + coefficients: dict + List of coefficients. If a coefficient is omitted, it is set to 0. + The valid coefficients are:\n + cL: str, callable, optional + Lift coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n + cQ: str, callable, optional + Side force coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n + cD: str, callable, optional + Drag coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n + cm: str, callable, optional + Pitch moment coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n + cn: str, callable, optional + Yaw moment coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n + cl: str, callable, optional + Roll moment coefficient. Can be a path to a CSV file or a callable. + Default is 0.\n center_of_pressure : tuple, list, optional Application point of the aerodynamic forces and moments. The center of pressure is defined in the local coordinate system of the @@ -76,6 +75,7 @@ def __init__( name : str, optional Name of the aerodynamic surface. Default is 'GenericSurface'. """ + self.reference_area = reference_area self.reference_length = reference_length self.center_of_pressure = center_of_pressure @@ -85,12 +85,80 @@ def __init__( self.cpz = center_of_pressure[2] self.name = name - self.cL = self._process_input(cL, "cL") - self.cD = self._process_input(cD, "cD") - self.cQ = self._process_input(cQ, "cQ") - self.cm = self._process_input(cm, "cm") - self.cn = self._process_input(cn, "cn") - self.cl = self._process_input(cl, "cl") + default_coefficients = self._get_default_coefficients() + self._check_coefficients(coefficients, default_coefficients) + coefficients = self._complete_coefficients(coefficients, default_coefficients) + for coeff, coeff_value in coefficients.items(): + value = self._process_input(coeff_value, coeff) + setattr(self, coeff, value) + + def _get_default_coefficients(self): + """Returns default coefficients + + Returns + ------- + default_coefficients: dict + Dictionary whose keys are the coefficients names and keys + are the default values. + """ + default_coefficients = { + "cL": 0, + "cQ": 0, + "cD": 0, + "cm": 0, + "cn": 0, + "cl": 0, + } + return default_coefficients + + def _complete_coefficients(self, input_coefficients, default_coefficients): + """Creates a copy of the input coefficients dict and fill it with missing + keys with default values + + Parameters + ---------- + input_coefficients : str, dict + Coefficients dictionary passed by the user. If the user only specifies some + of the coefficients, the remaining are completed with class default + values + default_coefficients : dict + Default coefficients of the class + + Returns + ------- + coefficients : dict + Coefficients dictionary used to setup coefficient attributes + """ + coefficients = copy.deepcopy(input_coefficients) + for coeff, value in default_coefficients.items(): + if coeff not in coefficients.keys(): + coefficients[coeff] = value + + return coefficients + + def _check_coefficients(self, input_coefficients, default_coefficients): + """Check if input coefficients have only valid keys + + Parameters + ---------- + input_coefficients : str, dict + Coefficients dictionary passed by the user. If the user only specifies some + of the coefficients, the remaining are completed with class default + values + default_coefficients : dict + Default coefficients of the class + + Raises + ------ + ValueError + Raises a value error if the input coefficient has an invalid key + """ + invalid_keys = set(input_coefficients) - set(default_coefficients) + if invalid_keys: + raise ValueError( + f"Invalid coefficient name(s) used in key(s): {', '.join(invalid_keys)}. " + "Check the documentation for valid names." + ) def _compute_from_coefficients( self, @@ -169,12 +237,8 @@ def compute_forces_and_moments( stream_mach, rho, cp, + omega, reynolds, - omega1, - omega2, - omega3, - *args, - **kwargs, ): """Computes the forces and moments acting on the aerodynamic surface. Used in each time step of the simulation. This method is valid for @@ -192,10 +256,10 @@ def compute_forces_and_moments( Air density. cp : Vector Center of pressure coordinates in the body frame. + omega: tuple[float, float, float] + Tuple containing angular velocities around the x, y, z axes. reynolds : float Reynolds number. - omega1, omega2, omega3 : float - Angular velocities around the x, y, z axes. Returns ------- @@ -218,9 +282,9 @@ def compute_forces_and_moments( beta, stream_mach, reynolds, - omega1, - omega2, - omega3, + omega[0], + omega[1], + omega[2], ) # Conversion from aerodynamic frame to body frame @@ -320,7 +384,7 @@ def __load_csv(self, file_path, coeff_name): reader = csv.reader(file) header = next(reader) except (FileNotFoundError, IOError) as e: - raise ValueError(f"Error reading {coeff_name} CSV file: {e}") + raise ValueError(f"Error reading {coeff_name} CSV file: {e}") from e if not header: raise ValueError(f"Invalid or empty CSV file for {coeff_name}.") diff --git a/rocketpy/rocket/aero_surface/linear_generic_surface.py b/rocketpy/rocket/aero_surface/linear_generic_surface.py index 863f999d9..35d5ae31e 100644 --- a/rocketpy/rocket/aero_surface/linear_generic_surface.py +++ b/rocketpy/rocket/aero_surface/linear_generic_surface.py @@ -11,42 +11,7 @@ def __init__( self, reference_area, reference_length, - cL_0=0, - cL_alpha=0, - cL_beta=0, - cL_p=0, - cL_q=0, - cL_r=0, - cQ_0=0, - cQ_alpha=0, - cQ_beta=0, - cQ_p=0, - cQ_q=0, - cQ_r=0, - cD_0=0, - cD_alpha=0, - cD_beta=0, - cD_p=0, - cD_q=0, - cD_r=0, - cm_0=0, - cm_alpha=0, - cm_beta=0, - cm_p=0, - cm_q=0, - cm_r=0, - cn_0=0, - cn_alpha=0, - cn_beta=0, - cn_p=0, - cn_q=0, - cn_r=0, - cl_0=0, - cl_alpha=0, - cl_beta=0, - cl_p=0, - cl_q=0, - cl_r=0, + coefficients, center_of_pressure=(0, 0, 0), name="Generic Linear Surface", ): @@ -75,112 +40,115 @@ def __init__( reference_length : int, float Reference length of the aerodynamic surface. Has the unit of meters. Commonly defined as the rocket's diameter. - cL_0 : callable, str, optional - Coefficient of lift at zero angle of attack. Default is 0. - cL_alpha : callable, str, optional - Coefficient of lift derivative with respect to angle of attack. - Default is 0. - cL_beta : callable, str, optional - Coefficient of lift derivative with respect to sideslip angle. - Default is 0. - cL_p : callable, str, optional - Coefficient of lift derivative with respect to roll rate. - Default is 0. - cL_q : callable, str, optional - Coefficient of lift derivative with respect to pitch rate. - Default is 0. - cL_r : callable, str, optional - Coefficient of lift derivative with respect to yaw rate. - Default is 0. - cQ_0 : callable, str, optional - Coefficient of pitch moment at zero angle of attack. - Default is 0. - cQ_alpha : callable, str, optional - Coefficient of pitch moment derivative with respect to angle of - attack. Default is 0. - cQ_beta : callable, str, optional - Coefficient of pitch moment derivative with respect to sideslip - angle. Default is 0. - cQ_p : callable, str, optional - Coefficient of pitch moment derivative with respect to roll rate. - Default is 0. - cQ_q : callable, str, optional - Coefficient of pitch moment derivative with respect to pitch rate. - Default is 0. - cQ_r : callable, str, optional - Coefficient of pitch moment derivative with respect to yaw rate. - Default is 0. - cD_0 : callable, str, optional - Coefficient of drag at zero angle of attack. Default is 0. - cD_alpha : callable, str, optional - Coefficient of drag derivative with respect to angle of attack. - Default is 0. - cD_beta : callable, str, optional - Coefficient of drag derivative with respect to sideslip angle. - Default is 0. - cD_p : callable, str, optional - Coefficient of drag derivative with respect to roll rate. - Default is 0. - cD_q : callable, str, optional - Coefficient of drag derivative with respect to pitch rate. - Default is 0. - cD_r : callable, str, optional - Coefficient of drag derivative with respect to yaw rate. - Default is 0. - cm_0 : callable, str, optional - Coefficient of pitch moment at zero angle of attack. - Default is 0. - cm_alpha : callable, str, optional - Coefficient of pitch moment derivative with respect to angle of - attack. Default is 0. - cm_beta : callable, str, optional - Coefficient of pitch moment derivative with respect to sideslip - angle. Default is 0. - cm_p : callable, str, optional - Coefficient of pitch moment derivative with respect to roll rate. - Default is 0. - cm_q : callable, str, optional - Coefficient of pitch moment derivative with respect to pitch rate. - Default is 0. - cm_r : callable, str, optional - Coefficient of pitch moment derivative with respect to yaw rate. - Default is 0. - cn_0 : callable, str, optional - Coefficient of yaw moment at zero angle of attack. - Default is 0. - cn_alpha : callable, str, optional - Coefficient of yaw moment derivative with respect to angle of - attack. Default is 0. - cn_beta : callable, str, optional - Coefficient of yaw moment derivative with respect to sideslip angle. - Default is 0. - cn_p : callable, str, optional - Coefficient of yaw moment derivative with respect to roll rate. - Default is 0. - cn_q : callable, str, optional - Coefficient of yaw moment derivative with respect to pitch rate. - Default is 0. - cn_r : callable, str, optional - Coefficient of yaw moment derivative with respect to yaw rate. - Default is 0. - cl_0 : callable, str, optional - Coefficient of roll moment at zero angle of attack. - Default is 0. - cl_alpha : callable, str, optional - Coefficient of roll moment derivative with respect to angle of - attack. Default is 0. - cl_beta : callable, str, optional - Coefficient of roll moment derivative with respect to sideslip - angle. Default is 0. - cl_p : callable, str, optional - Coefficient of roll moment derivative with respect to roll rate. - Default is 0. - cl_q : callable, str, optional - Coefficient of roll moment derivative with respect to pitch rate. - Default is 0. - cl_r : callable, str, optional - Coefficient of roll moment derivative with respect to yaw rate. - Default is 0. + coefficients: dict, optional + List of coefficients. If a coefficient is omitted, it is set to 0. + The valid coefficients are:\n + cL_0: callable, str, optional + Coefficient of lift at zero angle of attack. Default is 0.\n + cL_alpha: callable, str, optional + Coefficient of lift derivative with respect to angle of attack. + Default is 0.\n + cL_beta: callable, str, optional + Coefficient of lift derivative with respect to sideslip angle. + Default is 0.\n + cL_p: callable, str, optional + Coefficient of lift derivative with respect to roll rate. + Default is 0.\n + cL_q: callable, str, optional + Coefficient of lift derivative with respect to pitch rate. + Default is 0.\n + cL_r: callable, str, optional + Coefficient of lift derivative with respect to yaw rate. + Default is 0.\n + cQ_0: callable, str, optional + Coefficient of pitch moment at zero angle of attack. + Default is 0.\n + cQ_alpha: callable, str, optional + Coefficient of pitch moment derivative with respect to angle of + attack. Default is 0.\n + cQ_beta: callable, str, optional + Coefficient of pitch moment derivative with respect to sideslip + angle. Default is 0.\n + cQ_p: callable, str, optional + Coefficient of pitch moment derivative with respect to roll rate. + Default is 0.\n + cQ_q: callable, str, optional + Coefficient of pitch moment derivative with respect to pitch rate. + Default is 0.\n + cQ_r: callable, str, optional + Coefficient of pitch moment derivative with respect to yaw rate. + Default is 0.\n + cD_0: callable, str, optional + Coefficient of drag at zero angle of attack. Default is 0.\n + cD_alpha: callable, str, optional + Coefficient of drag derivative with respect to angle of attack. + Default is 0.\n + cD_beta: callable, str, optional + Coefficient of drag derivative with respect to sideslip angle. + Default is 0.\n + cD_p: callable, str, optional + Coefficient of drag derivative with respect to roll rate. + Default is 0.\n + cD_q: callable, str, optional + Coefficient of drag derivative with respect to pitch rate. + Default is 0.\n + cD_r: callable, str, optional + Coefficient of drag derivative with respect to yaw rate. + Default is 0.\n + cm_0: callable, str, optional + Coefficient of pitch moment at zero angle of attack. + Default is 0.\n + cm_alpha: callable, str, optional + Coefficient of pitch moment derivative with respect to angle of + attack. Default is 0.\n + cm_beta: callable, str, optional + Coefficient of pitch moment derivative with respect to sideslip + angle. Default is 0.\n + cm_p: callable, str, optional + Coefficient of pitch moment derivative with respect to roll rate. + Default is 0.\n + cm_q: callable, str, optional + Coefficient of pitch moment derivative with respect to pitch rate. + Default is 0.\n + cm_r: callable, str, optional + Coefficient of pitch moment derivative with respect to yaw rate. + Default is 0.\n + cn_0: callable, str, optional + Coefficient of yaw moment at zero angle of attack. + Default is 0.\n + cn_alpha: callable, str, optional + Coefficient of yaw moment derivative with respect to angle of + attack. Default is 0.\n + cn_beta: callable, str, optional + Coefficient of yaw moment derivative with respect to sideslip angle. + Default is 0.\n + cn_p: callable, str, optional + Coefficient of yaw moment derivative with respect to roll rate. + Default is 0.\n + cn_q: callable, str, optional + Coefficient of yaw moment derivative with respect to pitch rate. + Default is 0.\n + cn_r: callable, str, optional + Coefficient of yaw moment derivative with respect to yaw rate. + Default is 0.\n + cl_0: callable, str, optional + Coefficient of roll moment at zero angle of attack. + Default is 0.\n + cl_alpha: callable, str, optional + Coefficient of roll moment derivative with respect to angle of + attack. Default is 0.\n + cl_beta: callable, str, optional + Coefficient of roll moment derivative with respect to sideslip + angle. Default is 0.\n + cl_p: callable, str, optional + Coefficient of roll moment derivative with respect to roll rate. + Default is 0.\n + cl_q: callable, str, optional + Coefficient of roll moment derivative with respect to pitch rate. + Default is 0.\n + cl_r: callable, str, optional + Coefficient of roll moment derivative with respect to yaw rate. + Default is 0.\n center_of_pressure : tuple, optional Application point of the aerodynamic forces and moments. The center of pressure is defined in the local coordinate system of the @@ -188,54 +156,66 @@ def __init__( name : str Name of the aerodynamic surface. Default is 'GenericSurface'. """ - self.reference_area = reference_area - self.reference_length = reference_length - self.center_of_pressure = center_of_pressure - self.cp = center_of_pressure - self.cpx = center_of_pressure[0] - self.cpy = center_of_pressure[1] - self.cpz = center_of_pressure[2] - self.name = name - - self.cL_0 = self._process_input(cL_0, "cL_0") - self.cL_alpha = self._process_input(cL_alpha, "cL_alpha") - self.cL_beta = self._process_input(cL_beta, "cL_beta") - self.cL_p = self._process_input(cL_p, "cL_p") - self.cL_q = self._process_input(cL_q, "cL_q") - self.cL_r = self._process_input(cL_r, "cL_r") - self.cQ_0 = self._process_input(cQ_0, "cQ_0") - self.cQ_alpha = self._process_input(cQ_alpha, "cQ_alpha") - self.cQ_beta = self._process_input(cQ_beta, "cQ_beta") - self.cQ_p = self._process_input(cQ_p, "cQ_p") - self.cQ_q = self._process_input(cQ_q, "cQ_q") - self.cQ_r = self._process_input(cQ_r, "cQ_r") - self.cD_0 = self._process_input(cD_0, "cD_0") - self.cD_alpha = self._process_input(cD_alpha, "cD_alpha") - self.cD_beta = self._process_input(cD_beta, "cD_beta") - self.cD_p = self._process_input(cD_p, "cD_p") - self.cD_q = self._process_input(cD_q, "cD_q") - self.cD_r = self._process_input(cD_r, "cD_r") - self.cl_0 = self._process_input(cl_0, "cl_0") - self.cl_alpha = self._process_input(cl_alpha, "cl_alpha") - self.cl_beta = self._process_input(cl_beta, "cl_beta") - self.cl_p = self._process_input(cl_p, "cl_p") - self.cl_q = self._process_input(cl_q, "cl_q") - self.cl_r = self._process_input(cl_r, "cl_r") - self.cm_0 = self._process_input(cm_0, "cm_0") - self.cm_alpha = self._process_input(cm_alpha, "cm_alpha") - self.cm_beta = self._process_input(cm_beta, "cm_beta") - self.cm_p = self._process_input(cm_p, "cm_p") - self.cm_q = self._process_input(cm_q, "cm_q") - self.cm_r = self._process_input(cm_r, "cm_r") - self.cn_0 = self._process_input(cn_0, "cn_0") - self.cn_alpha = self._process_input(cn_alpha, "cn_alpha") - self.cn_beta = self._process_input(cn_beta, "cn_beta") - self.cn_p = self._process_input(cn_p, "cn_p") - self.cn_q = self._process_input(cn_q, "cn_q") - self.cn_r = self._process_input(cn_r, "cn_r") + + super().__init__( + reference_area=reference_area, + reference_length=reference_length, + coefficients=coefficients, + center_of_pressure=center_of_pressure, + name=name, + ) self.compute_all_coefficients() + def _get_default_coefficients(self): + """Returns default coefficients + + Returns + ------- + default_coefficients: dict + Dictionary whose keys are the coefficients names and keys + are the default values. + """ + default_coefficients = { + "cL_0": 0, + "cL_alpha": 0, + "cL_beta": 0, + "cL_p": 0, + "cL_q": 0, + "cL_r": 0, + "cQ_0": 0, + "cQ_alpha": 0, + "cQ_beta": 0, + "cQ_p": 0, + "cQ_q": 0, + "cQ_r": 0, + "cD_0": 0, + "cD_alpha": 0, + "cD_beta": 0, + "cD_p": 0, + "cD_q": 0, + "cD_r": 0, + "cm_0": 0, + "cm_alpha": 0, + "cm_beta": 0, + "cm_p": 0, + "cm_q": 0, + "cm_r": 0, + "cn_0": 0, + "cn_alpha": 0, + "cn_beta": 0, + "cn_p": 0, + "cn_q": 0, + "cn_r": 0, + "cl_0": 0, + "cl_alpha": 0, + "cl_beta": 0, + "cl_p": 0, + "cl_q": 0, + "cl_r": 0, + } + return default_coefficients + def compute_forcing_coefficient(self, c_0, c_alpha, c_beta): """Compute the forcing coefficient from the derivatives of the aerodynamic coefficients.""" @@ -297,6 +277,7 @@ def total_coefficient( def compute_all_coefficients(self): """Compute all the aerodynamic coefficients from the derivatives.""" + # pylint: disable=invalid-name self.cLf = self.compute_forcing_coefficient( self.cL_0, self.cL_alpha, self.cL_beta ) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 0eba7308e..6fafc1e85 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1464,16 +1464,18 @@ def u_dot( a32 = 2 * (e2 * e3 + e0 * e1) a33 = 1 - 2 * (e1**2 + e2**2) # Transformation matrix: (123) -> (XYZ) - K = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]] + K = Matrix([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]) + Kt = K.transpose # Calculate Forces and Moments # Get freestream speed wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) + speed_of_sound = self.env.speed_of_sound.get_value_opt(z) free_stream_speed = ( (wind_velocity_x - vx) ** 2 + (wind_velocity_y - vy) ** 2 + (vz) ** 2 ) ** 0.5 - free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z) + free_stream_mach = free_stream_speed / speed_of_sound # Determine aerodynamics forces # Determine Drag Force @@ -1507,76 +1509,47 @@ def u_dot( vy_b = a12 * vx + a22 * vy + a32 * vz vz_b = a13 * vx + a23 * vy + a33 * vz # Calculate lift and moment for each component of the rocket - for aero_surface, position in self.rocket.aerodynamic_surfaces: - comp_cp = ( - position - self.rocket.center_of_dry_mass_position - ) * self.rocket._csys - aero_surface.cpz - reference_area = aero_surface.reference_area - reference_length = aero_surface.reference_length + velocity_in_body_frame = Vector([vx_b, vy_b, vz_b]) + w = Vector([omega1, omega2, omega3]) + for aero_surface, _ in self.rocket.aerodynamic_surfaces: + # Component cp relative to CDM in body frame + comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface] # Component absolute velocity in body frame - comp_vx_b = vx_b + comp_cp * omega2 - comp_vy_b = vy_b - comp_cp * omega1 - comp_vz_b = vz_b - # Wind velocity at component - comp_z = z + comp_cp + comp_vb = velocity_in_body_frame + (w ^ comp_cp) + # Wind velocity at component altitude + comp_z = z + (K @ comp_cp).z comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z) comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z) # Component freestream velocity in body frame - comp_wind_vx_b = a11 * comp_wind_vx + a21 * comp_wind_vy - comp_wind_vy_b = a12 * comp_wind_vx + a22 * comp_wind_vy - comp_wind_vz_b = a13 * comp_wind_vx + a23 * comp_wind_vy - comp_stream_vx_b = comp_wind_vx_b - comp_vx_b - comp_stream_vy_b = comp_wind_vy_b - comp_vy_b - comp_stream_vz_b = comp_wind_vz_b - comp_vz_b - comp_stream_speed = ( - comp_stream_vx_b**2 + comp_stream_vy_b**2 + comp_stream_vz_b**2 - ) ** 0.5 - # Component attack angle and lift force - comp_attack_angle = 0 - comp_lift, comp_lift_xb, comp_lift_yb = 0, 0, 0 - if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0: - # normalize component stream velocity in body frame - comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed - if -1 * comp_stream_vz_bn < 1: - comp_attack_angle = np.arccos(-comp_stream_vz_bn) - c_lift = aero_surface.cl.get_value_opt( - comp_attack_angle, free_stream_mach - ) - # component lift force magnitude - comp_lift = ( - 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift - ) - # component lift force components - lift_dir_norm = (comp_stream_vx_b**2 + comp_stream_vy_b**2) ** 0.5 - comp_lift_xb = comp_lift * (comp_stream_vx_b / lift_dir_norm) - comp_lift_yb = comp_lift * (comp_stream_vy_b / lift_dir_norm) - # add to total lift force - R1 += comp_lift_xb - R2 += comp_lift_yb - # Add to total moment - M1 -= (comp_cp) * comp_lift_yb - M2 += (comp_cp) * comp_lift_xb - # Calculates Roll Moment - try: - clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters - M3_forcing = ( - (1 / 2 * rho * free_stream_speed**2) - * reference_area - * reference_length - * clf_delta.get_value_opt(free_stream_mach) - * cant_angle_rad - ) - M3_damping = ( - (1 / 2 * rho * free_stream_speed) - * reference_area - * (reference_length) ** 2 - * cld_omega.get_value_opt(free_stream_mach) - * omega3 - / 2 - ) - M3 += M3_forcing - M3_damping - except AttributeError: - pass + comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) + comp_stream_velocity = comp_wind_vb - comp_vb + comp_stream_speed = abs(comp_stream_velocity) + comp_stream_mach = comp_stream_speed / speed_of_sound + # Reynolds at component altitude + # TODO: Reynolds is only used in generic surfaces. This calculation + # should be moved to the surface class for efficiency + comp_reynolds = ( + self.env.density.get_value_opt(comp_z) + * comp_stream_speed + * aero_surface.reference_length + / self.env.dynamic_viscosity.get_value_opt(comp_z) + ) + # Forces and moments + X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( + comp_stream_velocity, + comp_stream_speed, + comp_stream_mach, + rho, + comp_cp, + w, + comp_reynolds, + ) + R1 += X + R2 += Y + R3 += Z + M1 += M + M2 += N + M3 += L # Off center moment M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1 @@ -1663,7 +1636,7 @@ def u_dot( (R3 - b * propellant_mass_at_t * (alpha2 - omega1 * omega3) + thrust) / total_mass_at_t, ] - ax, ay, az = np.dot(K, L) + ax, ay, az = K @ Vector(L) az -= self.env.gravity.get_value_opt(z) # Include gravity # Create u_dot @@ -1806,10 +1779,10 @@ def u_dot_generalized( # TODO: Reynolds is only used in generic surfaces. This calculation # should be moved to the surface class for efficiency comp_reynolds = ( - self.env.density.get_value_opt(z) + self.env.density.get_value_opt(comp_z) * comp_stream_speed * aero_surface.reference_length - / self.env.dynamic_viscosity.get_value_opt(z) + / self.env.dynamic_viscosity.get_value_opt(comp_z) ) # Forces and moments X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( @@ -1818,10 +1791,8 @@ def u_dot_generalized( comp_stream_mach, rho, comp_cp, + w, comp_reynolds, - omega1, - omega2, - omega3, ) R1 += X R2 += Y