diff --git a/.pylintrc b/.pylintrc index 2dddcb3bf..787f73e51 100644 --- a/.pylintrc +++ b/.pylintrc @@ -359,6 +359,9 @@ max-statements=25 # Minimum number of public methods for a class (see R0903). min-public-methods=0 +# Maximum number of positional arguments (see R0917). +max-positional-arguments = 25 + [EXCEPTIONS] diff --git a/CHANGELOG.md b/CHANGELOG.md index 49b65db25..3e5139933 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,11 +32,13 @@ Attention: The newest changes should be on top --> ### Added -- +- ENH: Generic Surfaces and Generic Linear Surfaces [#680](https://github.com/RocketPy-Team/RocketPy/pull/680) +- ENH: Free-Form Fins [#694](https://github.com/RocketPy-Team/RocketPy/pull/694) +- ENH: Expand Polation Options for ND Functions. [#691](https://github.com/RocketPy-Team/RocketPy/pull/691) ### Changed -- ENH: Expand Polation Options for ND Functions. [#691](https://github.com/RocketPy-Team/RocketPy/pull/691) +- ### Fixed diff --git a/README.md b/README.md index 5808678d3..c549840f7 100644 --- a/README.md +++ b/README.md @@ -329,6 +329,9 @@ RocketPy extends its gratitude to the following institutions for their support a Space Enterprise at Berkeley Logo + + Faraday Rocketry UPV Logo + ## Individual Contributors diff --git a/docs/reference/classes/aero_surfaces/GenericSurface.rst b/docs/reference/classes/aero_surfaces/GenericSurface.rst new file mode 100644 index 000000000..ed57c23af --- /dev/null +++ b/docs/reference/classes/aero_surfaces/GenericSurface.rst @@ -0,0 +1,5 @@ +Generic Surface Class +--------------------- + +.. autoclass:: rocketpy.GenericSurface + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/LinearGenericSurface.rst b/docs/reference/classes/aero_surfaces/LinearGenericSurface.rst new file mode 100644 index 000000000..0c39469fb --- /dev/null +++ b/docs/reference/classes/aero_surfaces/LinearGenericSurface.rst @@ -0,0 +1,5 @@ +Linear Generic Surface Class +---------------------------- + +.. autoclass:: rocketpy.LinearGenericSurface + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/index.rst b/docs/reference/classes/aero_surfaces/index.rst index 5cbb9eedb..c6e6a3efe 100644 --- a/docs/reference/classes/aero_surfaces/index.rst +++ b/docs/reference/classes/aero_surfaces/index.rst @@ -13,3 +13,5 @@ AeroSurface Classes EllipticalFins RailButtons AirBrakes + GenericSurface + LinearGenericSurface diff --git a/docs/static/institutional/faraday_team_logo.jpg b/docs/static/institutional/faraday_team_logo.jpg new file mode 100644 index 000000000..1ef1c2b85 Binary files /dev/null and b/docs/static/institutional/faraday_team_logo.jpg differ diff --git a/docs/static/rocket/aeroframe.png b/docs/static/rocket/aeroframe.png new file mode 100644 index 000000000..c48c86e1a Binary files /dev/null and b/docs/static/rocket/aeroframe.png differ diff --git a/docs/user/rocket/generic_surface.rst b/docs/user/rocket/generic_surface.rst new file mode 100644 index 000000000..8b0784a42 --- /dev/null +++ b/docs/user/rocket/generic_surface.rst @@ -0,0 +1,436 @@ +.. _genericsurfaces: + +Generic Surfaces and Custom Aerodynamic Coefficients +==================================================== + +Generic aerodynamic surfaces can be used to model aerodynamic forces based on +force and moment coefficients. The :class:`rocketpy.GenericSurface` receives the +coefficients as functions of the angle of attack, side slip angle, Mach number, +Reynolds number, pitch rate, yaw rate, and roll rate. + +The :class:`rocketpy.LinearGenericSurface` class model aerodynamic forces based +the force and moment coefficients derivatives. The coefficients are derivatives +of the force and moment coefficients with respect to the angle of attack, side +slip angle, Mach number, Reynolds number, pitch rate, yaw rate, and roll rate. + +These classes allows the user to be less dependent on the built-in aerodynamic +surfaces and to define their own aerodynamic coefficients. + +Both classes base their coefficient on the definition of the aerodynamic frame +of reference. + +Aerodynamic Frame +----------------- + +The aerodynamic frame of reference of the rocket is defined as follows: + +- The origin is at the rocket's center of dry mass (``center_of_dry_mass_position``). +- The ``z`` axis is defined along the rocket's centerline, pointing from the center of dry mass towards the nose. +- The ``x`` and ``y`` axes are perpendicular. +- The partial angle of attack (``alpha``) is defined as the angle, in the y-z + plane, from the velocity vector to the z axis. +- The partial side slip angle (``beta``) is defined as the angle, in the x-z + plane, from the velocity vector to the z axis. + +The following figure shows the aerodynamic frame of reference: + +.. figure:: ../../static/rocket/aeroframe.png + :align: center + :alt: Aerodynamic frame of reference + +In the figure we define: + +- :math:`\mathbf{\vec{V}}` as rocket velocity vector. +- :math:`x_B`, :math:`y_B`, and :math:`z_B` as the body axes. +- :math:`x_A`, :math:`y_A`, and :math:`z_A` as the aerodynamic axes. +- :math:`\alpha` as the partial angle of attack. +- :math:`\beta` as the side slip angle. +- :math:`L` as the lift force. +- :math:`D` as the drag force. +- :math:`Q` as the side force. + +Here we define the aerodynamic forces in the aerodynamic frame of reference as: + +.. math:: + \vec{\mathbf{F}}_A=\begin{bmatrix}X_A\\Y_A\\Z_A\end{bmatrix}_A=\begin{bmatrix}Q\\-L\\-D\end{bmatrix}_A + +The aerodynamic forces in the body axes coordinate system are defined as +:math:`\vec{\mathbf{F}}_B`. + +.. math:: + \vec{\mathbf{F}}_B=\begin{bmatrix}X_A\\Y_A\\Z_A\end{bmatrix}_B=\mathbf{M}_{BA}\cdot\begin{bmatrix}Q\\-L\\-D\end{bmatrix}_A + +Where the transformation matrix :math:`\mathbf{M}_{BA}`, which transforms the +aerodynamic forces from the aerodynamic frame of reference to the body axes +coordinate system, is defined as: + +.. math:: + \mathbf{M}_{BA} = \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & \cos(\alpha) & -\sin(\alpha) \\ + 0 & \sin(\alpha) & \cos(\alpha) + \end{bmatrix} + \begin{bmatrix} + \cos(\beta) & 0 & -\sin(\beta) \\ + 0 & 1 & 0 \\ + \sin(\beta) & 0 & \cos(\beta) + \end{bmatrix} + + +The forces coefficients can finally be defined as: + +- :math:`C_L` as the lift coefficient. +- :math:`C_Q` as the side force coefficient (or cross stream force coefficient). +- :math:`C_D` as the drag coefficient. + +And the forces from the coefficients are defined as: + +.. math:: + \begin{bmatrix}X_A\\Y_A\\Z_A\end{bmatrix}_B =\mathbf{M}_{BA}\cdot\overline{q}\cdot A_{ref}\cdot\begin{bmatrix}C_Q\\-C_L\\-C_D\end{bmatrix}_A + +Where: + +- :math:`\bar{q}` is the dynamic pressure. +- :math:`A_{ref}` is the reference area used to calculate the coefficients. + Commonly the rocket's cross-sectional area is used as the reference area. + +The moment coefficients can be defined as: + +- :math:`C_l` as the rolling moment coefficient. +- :math:`C_m` as the pitching moment coefficient. +- :math:`C_n` as the yawing moment coefficient. + +And the moments from the coefficients are defined as: + +.. math:: + \vec{\mathbf{M}}_B=\begin{bmatrix}M_{x_A}\\M_{y_A}\\M_{z_A}\end{bmatrix}_B =\overline{q}\cdot A_{ref}\cdot L_{ref}\cdot\begin{bmatrix}C_m\\C_n\\C_l\end{bmatrix} + +Where: + +- :math:`L_{ref}` is the reference length used to calculate the coefficients. + Commonly the rocket's diameter is used as the reference length. + + +Aerodynamic angles +~~~~~~~~~~~~~~~~~~ + +There aerodynamic angles are defined in two different ways in RocketPy: + +- As the angle of attack (:math:`\alpha`) and the side slip \ + angle (:math:`\beta`), which are defined in the image above. These are used \ + in the calculation of the generic surface forces and moments. +- As the total angle of attack (:math:`\alpha_{\text{tot}}`), defined as the \ + angle between the total velocity vector and the rocket's centerline. This is \ + used in the calculation of the standard aerodynamic surface forces and moments. + +The partial angles are calculated as: + +.. math:: + \begin{aligned} + \alpha &= \arctan\left(\frac{V_y}{V_z}\right) \\ + \beta &= \arctan\left(\frac{V_x}{V_z}\right) + \end{aligned} + +The total angle of attack is calculated as: + +.. math:: + \alpha_{\text{tot}} = \arccos\left(\frac{\mathbf{\vec{V}}\cdot\mathbf{z_B}}{||\mathbf{\vec{V}}||\cdot||\mathbf{z_B}||}\right) + +.. note:: + When the simulation is done, the total angle of attack is accessed through + the :attr:`rocketpy.Flight.angle_of_attack` attribute. + The partial angles of attack and side slip are accessed through the + :attr:`rocketpy.Flight.partial_angle_of_attack` and + :attr:`rocketpy.Flight.angle_of_sideslip` attributes, respectively. + +.. _genericsurface: + +Generic Surface Class +--------------------- + +The :class:`rocketpy.GenericSurface` class is used to define an aerodynamic +surface based on force and moment coefficients. A generic surface is defined +as follows: + +.. seealso:: + For more information on class initialization, see + :class:`rocketpy.GenericSurface.__init__` + + +.. code-block:: python + + from rocketpy import GenericSurface + + radius = 0.0635 + + generic_surface = GenericSurface( + reference_area=np.pi * radius**2, + reference_length=2 * radius, + coefficients={ + "cL": "cL.csv", + "cQ": "cQ.csv", + "cD": "cD.csv", + "cm": "cm.csv", + "cn": "cn.csv", + "cl": "cl.csv", + }, + name="Generic Surface", + ) + +The ``coefficients`` argument is a dictionary containing the coefficients of the +generic surface. The keys of the dictionary are the coefficient names, and the +values are the coefficients. The possible coefficient names are: + +- ``cL``: Lift coefficient. +- ``cQ``: Side force coefficient. +- ``cD``: Drag coefficient. +- ``cm``: Pitching moment coefficient. +- ``cn``: Yawing moment coefficient. +- ``cl``: Rolling moment coefficient. + +Only one of the coefficients is required to be provided, but any combination of +the coefficients can be used. The coefficient values can be provided as a +single value, a callable function of seven arguments, or a path to a ``.csv`` +file containing the values. + +The coefficients are all functions of: + +- Angle of attack (:math:`\alpha`) in radians. +- Side slip angle (:math:`\beta`) in radians. +- Mach number (:math:`Ma`). +- Reynolds number (:math:`Re`). +- Pitch rate (:math:`q`) in radians per second. +- Yaw rate (:math:`r`) in radians per second. +- Roll rate (:math:`p`) in radians per second. + +.. math:: + \begin{aligned} + C_L &= f(\alpha, \beta, Ma, Re, q, r, p) \\ + C_Q &= f(\alpha, \beta, Ma, Re, q, r, p) \\ + C_D &= f(\alpha, \beta, Ma, Re, q, r, p) \\ + C_m &= f(\alpha, \beta, Ma, Re, q, r, p) \\ + C_n &= f(\alpha, \beta, Ma, Re, q, r, p) \\ + C_l &= f(\alpha, \beta, Ma, Re, q, r, p) + \end{aligned} + +From the coefficients, the forces and moments are calculated with + +.. math:: + \begin{aligned} + L &= \overline{q}\cdot A_{ref}\cdot C_L \\ + Q &= \overline{q}\cdot A_{ref}\cdot C_Q \\ + D &= \overline{q}\cdot A_{ref}\cdot C_D \\ + M_{m} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_m \\ + M_{n} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_n \\ + M_{l} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_l + \end{aligned} + +These coefficients can be defined as a callable such as: + +.. code-block:: python + + def coefficient(alpha, beta, Ma, Re, q, r, p): + ... + return value + +In which any algorithm can be implemented to calculate the coefficient values. + +Otherwise, the coefficients can be defined as a ``.csv`` file. The file must +contain a header with at least one of the following columns representing the +independent variables: + +- ``alpha``: Angle of attack. +- ``beta``: Side slip angle. +- ``mach``: Mach number. +- ``reynolds``: Reynolds number. +- ``q``: Pitch rate. +- ``r``: Yaw rate. +- ``p``: Roll rate. + +The last column must be the coefficient value, and must contain a header, +though the header name can be anything. + +.. important:: + Not all columns need to be present in the file, but the columns that are + present must be named, **and ordered**, as described above. + +An example of a ``.csv`` file is shown below: + +.. code-block:: + + "alpha", "mach", "coefficient" + -0.017, 0, -0.11 + -0.017, 1, -0.127 + -0.017, 2, -0.084 + -0.017, 3, -0.061 + 0.0, 0, 0.0 + 0.0, 1, 0.0 + 0.0, 2, 0.0 + 0.0, 3, 0.0 + 0.017, 0, 0.11 + 0.017, 1, 0.127 + 0.017, 2, 0.084 + 0.017, 3, 0.061 + +After the definition of the ``GenericSurface`` object, it must be added to the +rocket's configuration: + +.. seealso:: + For more information on how to add a generic surface to the rocket, see + :class:`rocketpy.Rocket.add_generic_surface` + +.. code-block:: python + :emphasize-lines: 5 + + from rocketpy import Rocket + rocket = Rocket( + ... + ) + rocket.add_surfaces(generic_surface, position=(0,0,0)) + +The position of the generic surface is defined in the User Defined coordinate +System, see :ref:`rocketaxes` for more information. + +.. tip:: + If defining the coefficients of the entire rocket is desired, only a single + generic surface can be added to the rocket, positioned at the center of dry + mass. This will be equivalent to defining the coefficients of the entire + rocket. + +.. attention:: + If there generic surface is positioned **not** at the center of dry mass, the + forces generated by the force coefficients (cL, cQ, cD) will generate a + moment around the center of dry mass. This moment will be calculated and + added to the moment generated by the moment coefficients (cm, cn, cl). + + +.. _lineargenericsurface: + +Linear Generic Surface Class +---------------------------- + +The :class:`rocketpy.LinearGenericSurface` class is used to define a aerodynamic +surface based on the forces and moments coefficient derivatives. A linear generic +surface will receive the derivatives of each coefficient with respect to the +independent variables. The derivatives are defined as: + +- :math:`C_{\alpha}=\frac{dC}{d\alpha}`: Coefficient derivative with respect to angle of attack. +- :math:`C_{\beta}=\frac{dC}{d\beta}`: Coefficient derivative with respect to side slip angle. +- :math:`C_{Ma}=\frac{dC}{dMa}`: Coefficient derivative with respect to Mach number. +- :math:`C_{Re}=\frac{dC}{dRe}`: Coefficient derivative with respect to Reynolds number. +- :math:`C_{q}=\frac{dC}{dq}`: Coefficient derivative with respect to pitch rate. +- :math:`C_{r}=\frac{dC}{dr}`: Coefficient derivative with respect to yaw rate. +- :math:`C_{p}=\frac{dC}{dp}`: Coefficient derivative with respect to roll rate. + +A non derivative coefficient :math:`C_{0}` is also included. + +Each coefficient derivative is defined as a function of all the seven +independent variables. + +The coefficients are then grouped into **forcing** coefficients: + +.. math:: + \begin{aligned} + C_{Lf} &= C_{L0} + C_{L\alpha}\cdot\alpha + C_{L\beta}\cdot\beta + C_{LMa}\cdot Ma + C_{LRe}\cdot Re \\ + C_{Qf} &= C_{Q0} + C_{Q\alpha}\cdot\alpha + C_{Q\beta}\cdot\beta + C_{QMa}\cdot Ma + C_{QRe}\cdot Re \\ + C_{Df} &= C_{D0} + C_{D\alpha}\cdot\alpha + C_{D\beta}\cdot\beta + C_{DMa}\cdot Ma + C_{DRe}\cdot Re \\ + C_{mf} &= C_{m0} + C_{m\alpha}\cdot\alpha + C_{m\beta}\cdot\beta + C_{mMa}\cdot Ma + C_{mRe}\cdot Re \\ + C_{nf} &= C_{n0} + C_{n\alpha}\cdot\alpha + C_{n\beta}\cdot\beta + C_{nMa}\cdot Ma + C_{nRe}\cdot Re \\ + C_{lf} &= C_{l0} + C_{l\alpha}\cdot\alpha + C_{l\beta}\cdot\beta + C_{lMa}\cdot Ma + C_{lRe}\cdot Re + \end{aligned} + +And **damping** coefficients: + +.. math:: + \begin{aligned} + C_{Ld} &= C_{L_{q}}\cdot q + C_{L_{r}}\cdot r + C_{L_{p}}\cdot p \\ + C_{Qd} &= C_{Q_{q}}\cdot q + C_{Q_{r}}\cdot r + C_{Q_{p}}\cdot p \\ + C_{Dd} &= C_{D_{q}}\cdot q + C_{D_{r}}\cdot r + C_{D_{p}}\cdot p \\ + C_{md} &= C_{m_{q}}\cdot q + C_{m_{r}}\cdot r + C_{m_{p}}\cdot p \\ + C_{nd} &= C_{n_{q}}\cdot q + C_{n_{r}}\cdot r + C_{n_{p}}\cdot p \\ + C_{ld} &= C_{l_{q}}\cdot q + C_{l_{r}}\cdot r + C_{l_{p}}\cdot p + \end{aligned} + +The forces and moments are then calculated as: + +.. math:: + \begin{aligned} + L &= \overline{q}\cdot A_{ref}\cdot C_{Lf} + \overline{q}\cdot A_{ref}\cdot \frac{L_{ref}}{2V} C_{Ld} \\ + Q &= \overline{q}\cdot A_{ref}\cdot C_{Qf} + \overline{q}\cdot A_{ref}\cdot \frac{L_{ref}}{2V} C_{Qd} \\ + D &= \overline{q}\cdot A_{ref}\cdot C_{Df} + \overline{q}\cdot A_{ref}\cdot \frac{L_{ref}}{2V} C_{Dd} \\ + M_{m} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_{mf} + \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot \frac{L_{ref}}{2V} C_{md} \\ + M_{n} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_{nf} + \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot \frac{L_{ref}}{2V} C_{nd} \\ + M_{l} &= \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot C_{lf} + \overline{q}\cdot A_{ref}\cdot L_{ref}\cdot \frac{L_{ref}}{2V} C_{ld} + \end{aligned} + +The linear generic surface is defined very similarly to the generic surface. +The coefficients are defined in the same way, but with the addition of the +derivative values. + +An example of a linear generic surface defined with **all** the coefficients is +shown below: + +.. seealso:: + For more information on class initialization, see + :class:`rocketpy.LinearGenericSurface.__init__` + +.. code-block:: python + + from rocketpy import LinearGenericSurface + linear_generic_surface = LinearGenericSurface( + reference_area=np.pi * 0.0635**2, + reference_length=2 * 0.0635, + coefficients={ + "cL_0": "cL_0.csv", + "cL_alpha": "cL_alpha.csv", + "cL_beta": "cL_beta.csv", + "cL_Ma": "cL_Ma.csv", + "cL_Re": "cL_Re.csv", + "cL_q": "cL_q.csv", + "cL_r": "cL_r.csv", + "cL_p": "cL_p.csv", + "cQ_0": "cQ_0.csv", + "cQ_alpha": "cQ_alpha.csv", + "cQ_beta": "cQ_beta.csv", + "cQ_Ma": "cQ_Ma.csv", + "cQ_Re": "cQ_Re.csv", + "cQ_q": "cQ_q.csv", + "cQ_r": "cQ_r.csv", + "cQ_p": "cQ_p.csv", + "cD_0": "cD_0.csv", + "cD_alpha": "cD_alpha.csv", + "cD_beta": "cD_beta.csv", + "cD_Ma": "cD_Ma.csv", + "cD_Re": "cD_Re.csv", + "cD_q": "cD_q.csv", + "cD_r": "cD_r.csv", + "cD_p": "cD_p.csv", + "cm_0": "cm_0.csv", + "cm_alpha": "cm_alpha.csv", + "cm_beta": "cm_beta.csv", + "cm_Ma": "cm_Ma.csv", + "cm_Re": "cm_Re.csv", + "cm_q": "cm_q.csv", + "cm_r": "cm_r.csv", + "cm_p": "cm_p.csv", + "cn_0": "cn_0.csv", + "cn_alpha": "cn_alpha.csv", + "cn_beta": "cn_beta.csv", + "cn_Ma": "cn_Ma.csv", + "cn_Re": "cn_Re.csv", + "cn_q": "cn_q.csv", + "cn_r": "cn_r.csv", + "cn_p": "cn_p.csv", + "cl_0": "cl_0.csv", + "cl_alpha": "cl_alpha.csv", + "cl_beta": "cl_beta.csv", + "cl_Ma": "cl_Ma.csv", + "cl_Re": "cl_Re.csv", + "cl_q": "cl_q.csv", + "cl_r": "cl_r.csv", + "cl_p": "cl_p.csv", + }, + ) + rocket.add_surfaces(linear_generic_surface, position=(0,0,0)) + + diff --git a/docs/user/rocket/rocket.rst b/docs/user/rocket/rocket.rst index 17ec63c1a..956175409 100644 --- a/docs/user/rocket/rocket.rst +++ b/docs/user/rocket/rocket.rst @@ -12,4 +12,10 @@ Rocket Usage :caption: Rocket Class Axes Definitions Rocket Class Axes Definitions + +.. toctree:: + :maxdepth: 3 + :caption: Generic Surfaces and Custom Aerodynamic Coefficients + + Generic Surfaces and Custom Aerodynamic Coefficients \ No newline at end of file diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 90ffcf72a..539b8b2cb 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -30,6 +30,9 @@ Components, EllipticalFins, Fins, + FreeFormFins, + GenericSurface, + LinearGenericSurface, NoseCone, Parachute, RailButtons, diff --git a/rocketpy/plots/aero_surface_plots.py b/rocketpy/plots/aero_surface_plots.py index c242973b3..2a696b1fc 100644 --- a/rocketpy/plots/aero_surface_plots.py +++ b/rocketpy/plots/aero_surface_plots.py @@ -380,6 +380,73 @@ def draw(self): plt.show() +class _FreeFormFinsPlots(_FinsPlots): + """Class that contains all free form fin plots.""" + + # pylint: disable=too-many-statements + 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 + """ + # Color cycle [#348ABD, #A60628, #7A68A6, #467821, #D55E00, #CC79A7, + # #56B4E9, #009E73, #F0E442, #0072B2] + + # Center of pressure + cp_point = [self.aero_surface.cpz, self.aero_surface.Yma] + + # Mean Aerodynamic Chord + yma_line = plt.Line2D( + ( + self.aero_surface.mac_lead, + self.aero_surface.mac_lead + self.aero_surface.mac_length, + ), + (self.aero_surface.Yma, self.aero_surface.Yma), + color="#467821", + linestyle="--", + label="Mean Aerodynamic Chord", + ) + + # Plotting + fig = plt.figure(figsize=(7, 4)) + with plt.style.context("bmh"): + ax = fig.add_subplot(111) + + # Fin + ax.scatter( + self.aero_surface.shape_vec[0], + self.aero_surface.shape_vec[1], + color="#A60628", + ) + ax.plot( + self.aero_surface.shape_vec[0], + self.aero_surface.shape_vec[1], + color="#A60628", + ) + # line from the last point to the first point + ax.plot( + [self.aero_surface.shape_vec[0][-1], self.aero_surface.shape_vec[0][0]], + [self.aero_surface.shape_vec[1][-1], self.aero_surface.shape_vec[1][0]], + color="#A60628", + ) + + ax.add_line(yma_line) + ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10) + ax.scatter(*cp_point, facecolors="none", edgecolors="red", s=500, zorder=10) + + # Plot settings + ax.set_xlabel("Root chord (m)") + ax.set_ylabel("Span (m)") + ax.set_title("Free Form Fin Cross Section") + ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + + plt.tight_layout() + plt.show() + + class _TailPlots(_AeroSurfacePlots): """Class that contains all tail plots.""" @@ -409,3 +476,17 @@ def all(self): None """ self.drag_coefficient_curve() + + +class _GenericSurfacePlots(_AeroSurfacePlots): + """Class that contains all generic surface plots.""" + + def draw(self): + pass + + +class _LinearGenericSurfacePlots(_AeroSurfacePlots): + """Class that contains all linear generic surface plots.""" + + def draw(self): + pass diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index 1e5c53449..0babaff4c 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -635,9 +635,9 @@ def fluid_mechanics_data(self): # pylint: disable=too-many-statements ------- None """ - plt.figure(figsize=(9, 12)) + plt.figure(figsize=(9, 16)) - ax1 = plt.subplot(411) + ax1 = plt.subplot(611) ax1.plot(self.flight.mach_number[:, 0], self.flight.mach_number[:, 1]) ax1.set_xlim(0, self.flight.t_final) ax1.set_title("Mach Number") @@ -645,7 +645,7 @@ def fluid_mechanics_data(self): # pylint: disable=too-many-statements ax1.set_ylabel("Mach Number") ax1.grid() - ax2 = plt.subplot(412) + ax2 = plt.subplot(612) ax2.plot(self.flight.reynolds_number[:, 0], self.flight.reynolds_number[:, 1]) ax2.set_xlim(0, self.flight.t_final) ax2.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) @@ -654,7 +654,7 @@ def fluid_mechanics_data(self): # pylint: disable=too-many-statements ax2.set_ylabel("Reynolds Number") ax2.grid() - ax3 = plt.subplot(413) + ax3 = plt.subplot(613) ax3.plot( self.flight.dynamic_pressure[:, 0], self.flight.dynamic_pressure[:, 1], @@ -678,7 +678,7 @@ def fluid_mechanics_data(self): # pylint: disable=too-many-statements ax3.set_ylabel("Pressure (Pa)") ax3.grid() - ax4 = plt.subplot(414) + ax4 = plt.subplot(614) ax4.plot(self.flight.angle_of_attack[:, 0], self.flight.angle_of_attack[:, 1]) ax4.set_title("Angle of Attack") ax4.set_xlabel("Time (s)") @@ -687,6 +687,33 @@ def fluid_mechanics_data(self): # pylint: disable=too-many-statements ax4.set_ylim(0, self.flight.angle_of_attack(self.flight.out_of_rail_time) + 15) ax4.grid() + ax5 = plt.subplot(615) + ax5.plot( + self.flight.partial_angle_of_attack[:, 0], + self.flight.partial_angle_of_attack[:, 1], + ) + ax5.set_title("Partial Angle of Attack") + ax5.set_xlabel("Time (s)") + ax5.set_ylabel("Partial Angle of Attack (°)") + ax5.set_xlim(self.flight.out_of_rail_time, self.first_event_time) + ax5.set_ylim( + 0, self.flight.partial_angle_of_attack(self.flight.out_of_rail_time) + 15 + ) + ax5.grid() + + ax6 = plt.subplot(616) + ax6.plot( + self.flight.angle_of_sideslip[:, 0], self.flight.angle_of_sideslip[:, 1] + ) + ax6.set_title("Angle of Sideslip") + ax6.set_xlabel("Time (s)") + ax6.set_ylabel("Angle of Sideslip (°)") + ax6.set_xlim(self.flight.out_of_rail_time, self.first_event_time) + ax6.set_ylim( + 0, self.flight.angle_of_sideslip(self.flight.out_of_rail_time) + 15 + ) + ax6.grid() + plt.subplots_adjust(hspace=0.5) plt.show() diff --git a/rocketpy/plots/monte_carlo_plots.py b/rocketpy/plots/monte_carlo_plots.py index 2264597da..7f8313d4b 100644 --- a/rocketpy/plots/monte_carlo_plots.py +++ b/rocketpy/plots/monte_carlo_plots.py @@ -54,9 +54,14 @@ def ellipses( "The image file was not found. Please check the path." ) from e - impact_ellipses, apogee_ellipses, apogee_x, apogee_y, impact_x, impact_y = ( - generate_monte_carlo_ellipses(self.monte_carlo.results) - ) + ( + impact_ellipses, + apogee_ellipses, + apogee_x, + apogee_y, + impact_x, + impact_y, + ) = generate_monte_carlo_ellipses(self.monte_carlo.results) # Create plot figure plt.figure(figsize=(8, 6), dpi=150) diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py index 2add1e78f..4ae47cdc4 100644 --- a/rocketpy/plots/rocket_plots.py +++ b/rocketpy/plots/rocket_plots.py @@ -5,6 +5,7 @@ from rocketpy.motors import EmptyMotor, HybridMotor, LiquidMotor, SolidMotor from rocketpy.rocket.aero_surface import Fins, NoseCone, Tail +from rocketpy.rocket.aero_surface.generic_surface import GenericSurface class _RocketPlots: @@ -208,7 +209,7 @@ def draw(self, vis_args=None, plane="xz"): reverse = csys == 1 self.rocket.aerodynamic_surfaces.sort_by_position(reverse=reverse) - drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args) + drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args, plane) last_radius, last_x = self._draw_tubes(ax, drawn_surfaces, vis_args) self._draw_motor(last_radius, last_x, ax, vis_args) self._draw_rail_buttons(ax, vis_args) @@ -224,7 +225,7 @@ def draw(self, vis_args=None, plane="xz"): plt.tight_layout() plt.show() - def _draw_aerodynamic_surfaces(self, ax, vis_args): + def _draw_aerodynamic_surfaces(self, ax, vis_args, plane): """Draws the aerodynamic surfaces and saves the position of the points of interest for the tubes.""" # List of drawn surfaces with the position of points of interest @@ -239,11 +240,15 @@ def _draw_aerodynamic_surfaces(self, ax, vis_args): for surface, position in self.rocket.aerodynamic_surfaces: if isinstance(surface, NoseCone): - self._draw_nose_cone(ax, surface, position, drawn_surfaces, vis_args) + self._draw_nose_cone(ax, surface, position.z, drawn_surfaces, vis_args) elif isinstance(surface, Tail): - self._draw_tail(ax, surface, position, drawn_surfaces, vis_args) + self._draw_tail(ax, surface, position.z, drawn_surfaces, vis_args) elif isinstance(surface, Fins): - self._draw_fins(ax, surface, position, drawn_surfaces, vis_args) + self._draw_fins(ax, surface, position.z, drawn_surfaces, vis_args) + elif isinstance(surface, GenericSurface): + self._draw_generic_surface( + ax, surface, position, drawn_surfaces, vis_args, plane + ) return drawn_surfaces def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args): @@ -336,6 +341,39 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): drawn_surfaces.append((surface, position, surface.rocket_radius, x_rotated[-1])) + def _draw_generic_surface( + self, + ax, + surface, + position, + drawn_surfaces, + vis_args, # pylint: disable=unused-argument + plane, + ): + """Draws the generic surface and saves the position of the points of interest + for the tubes.""" + if plane == "xz": + # z position of the sensor is the x position in the plot + x_pos = position[2] + # x position of the surface is the y position in the plot + y_pos = position[0] + elif plane == "yz": + # z position of the surface is the x position in the plot + x_pos = position[2] + # y position of the surface is the y position in the plot + y_pos = position[1] + else: + raise ValueError("Plane must be 'xz' or 'yz'.") + + ax.scatter( + x_pos, + y_pos, + linewidth=2, + zorder=10, + label=surface.name, + ) + drawn_surfaces.append((surface, position.z, self.rocket.radius, x_pos)) + def _draw_tubes(self, ax, drawn_surfaces, vis_args): """Draws the tubes between the aerodynamic surfaces.""" for i, d_surface in enumerate(drawn_surfaces): @@ -528,8 +566,8 @@ def _draw_rail_buttons(self, ax, vis_args): """Draws the rail buttons of the rocket.""" try: buttons, pos = self.rocket.rail_buttons[0] - lower = pos - upper = pos + buttons.buttons_distance * self.rocket._csys + lower = pos.z + upper = lower + buttons.buttons_distance * self.rocket._csys ax.scatter( lower, -self.rocket.radius, marker="s", color=vis_args["buttons"], s=15 ) diff --git a/rocketpy/plots/sensitivity_plots.py b/rocketpy/plots/sensitivity_plots.py index 1f507709c..023dd8f8e 100644 --- a/rocketpy/plots/sensitivity_plots.py +++ b/rocketpy/plots/sensitivity_plots.py @@ -3,7 +3,6 @@ class _SensitivityModelPlots: - def __init__(self, model): self.model = model diff --git a/rocketpy/prints/aero_surface_prints.py b/rocketpy/prints/aero_surface_prints.py index 7cc87c28f..378433511 100644 --- a/rocketpy/prints/aero_surface_prints.py +++ b/rocketpy/prints/aero_surface_prints.py @@ -72,7 +72,6 @@ def geometry(self): class _FinsPrints(_AeroSurfacePrints): - def geometry(self): print("Geometric information of the fin set:") print("-------------------------------------") @@ -184,6 +183,10 @@ class _EllipticalFinsPrints(_FinsPrints): """Class that contains all elliptical fins prints.""" +class _FreeFormFinsPrints(_FinsPrints): + """Class that contains all free form fins prints.""" + + class _TailPrints(_AeroSurfacePrints): """Class that contains all tail prints.""" @@ -228,3 +231,45 @@ def geometry(self): def all(self): pass + + +class _GenericSurfacePrints(_AeroSurfacePrints): + """Class that contains all generic surface prints.""" + + def geometry(self): + print("Geometric information of the Surface:") + print("----------------------------------") + print(f"Reference Area: {self.generic_surface.reference_area:.3f} m") + print(f"Reference length: {2*self.generic_surface.rocket_radius:.3f} m") + + def all(self): + """Prints all information of the generic surface. + + Returns + ------- + None + """ + self.identity() + self.geometry() + self.lift() + + +class _LinearGenericSurfacePrints(_AeroSurfacePrints): + """Class that contains all linear generic surface prints.""" + + def geometry(self): + print("Geometric information of the Surface:") + print("----------------------------------") + print(f"Reference Area: {self.generic_surface.reference_area:.3f} m") + print(f"Reference length: {2*self.generic_surface.rocket_radius:.3f} m") + + def all(self): + """Prints all information of the linear generic surface. + + Returns + ------- + None + """ + self.identity() + self.geometry() + self.lift() diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py index 7ad42fd7b..c9f5585ac 100644 --- a/rocketpy/prints/rocket_prints.py +++ b/rocketpy/prints/rocket_prints.py @@ -1,3 +1,6 @@ +from rocketpy.rocket.aero_surface.generic_surface import GenericSurface + + class _RocketPrints: """Class that holds prints methods for Rocket class. @@ -98,7 +101,9 @@ def rocket_aerodynamics_quantities(self): None """ print("\nAerodynamics Lift Coefficient Derivatives\n") - for surface, position in self.rocket.aerodynamic_surfaces: + for surface, _ in self.rocket.aerodynamic_surfaces: + if isinstance(surface, GenericSurface): + continue name = surface.name # ref_factor corrects lift for different reference areas ref_factor = (surface.rocket_radius / self.rocket.radius) ** 2 @@ -113,7 +118,7 @@ def rocket_aerodynamics_quantities(self): cpz = surface.cp[2] # relative to the user defined coordinate system print( f"{name} Center of Pressure position: " - f"{position - self.rocket._csys * cpz:.3f} m" + f"{position.z - self.rocket._csys * cpz:.3f} m" ) print("\nStability\n") print( diff --git a/rocketpy/prints/sensitivity_prints.py b/rocketpy/prints/sensitivity_prints.py index 014be586c..8a0301d79 100644 --- a/rocketpy/prints/sensitivity_prints.py +++ b/rocketpy/prints/sensitivity_prints.py @@ -4,7 +4,6 @@ class _SensitivityModelPrints: - def __init__(self, model): self.model = model diff --git a/rocketpy/prints/sensors_prints.py b/rocketpy/prints/sensors_prints.py index 73ab062f8..4bf844c5a 100644 --- a/rocketpy/prints/sensors_prints.py +++ b/rocketpy/prints/sensors_prints.py @@ -68,7 +68,6 @@ def all(self): class _InertialSensorPrints(_SensorPrints): - def orientation(self): """Prints the orientation of the sensor.""" print("\nOrientation of the Sensor:\n") diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py index fb6ea2b2c..0687b5ee5 100644 --- a/rocketpy/rocket/__init__.py +++ b/rocketpy/rocket/__init__.py @@ -4,6 +4,9 @@ AirBrakes, EllipticalFins, Fins, + FreeFormFins, + GenericSurface, + LinearGenericSurface, NoseCone, RailButtons, Tail, diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py index 9d9c68586..ad784f8d0 100644 --- a/rocketpy/rocket/aero_surface/__init__.py +++ b/rocketpy/rocket/aero_surface/__init__.py @@ -1,6 +1,13 @@ from rocketpy.rocket.aero_surface.aero_surface import AeroSurface from rocketpy.rocket.aero_surface.air_brakes import AirBrakes -from rocketpy.rocket.aero_surface.fins import EllipticalFins, Fins, TrapezoidalFins +from rocketpy.rocket.aero_surface.fins import ( + EllipticalFins, + Fins, + FreeFormFins, + TrapezoidalFins, +) +from rocketpy.rocket.aero_surface.generic_surface import GenericSurface +from rocketpy.rocket.aero_surface.linear_generic_surface import LinearGenericSurface from rocketpy.rocket.aero_surface.nose_cone import NoseCone from rocketpy.rocket.aero_surface.rail_buttons import RailButtons from rocketpy.rocket.aero_surface.tail import Tail diff --git a/rocketpy/rocket/aero_surface/aero_surface.py b/rocketpy/rocket/aero_surface/aero_surface.py index 81b5610e3..15ca14f1d 100644 --- a/rocketpy/rocket/aero_surface/aero_surface.py +++ b/rocketpy/rocket/aero_surface/aero_surface.py @@ -90,3 +90,60 @@ def all_info(self): ------- None """ + + def compute_forces_and_moments( + self, + stream_velocity, + stream_speed, + stream_mach, + rho, + cp, + *args, + ): # 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. + + Parameters + ---------- + stream_velocity : tuple + Tuple containing the stream velocity components in the body frame. + stream_speed : int, float + Speed of the stream in m/s. + stream_mach : int, float + Mach number of the stream. + rho : int, float + Density of the stream in kg/m^3. + cp : Vector + Center of pressure coordinates in the body frame. + args : tuple + Additional arguments. + kwargs : dict + Additional keyword arguments. + + 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, M3 = 0, 0, 0, 0, 0, 0 + cpz = cp[2] + stream_vx, stream_vy, stream_vz = stream_velocity + if stream_vx**2 + stream_vy**2 != 0: # TODO: maybe try/except + # Normalize component stream velocity in body frame + stream_vzn = stream_vz / stream_speed + if -1 * stream_vzn < 1: + attack_angle = np.arccos(-stream_vzn) + c_lift = self.cl.get_value_opt(attack_angle, stream_mach) + # Component lift force magnitude + lift = 0.5 * rho * (stream_speed**2) * self.reference_area * c_lift + # Component lift force components + lift_dir_norm = (stream_vx**2 + stream_vy**2) ** 0.5 + lift_xb = lift * (stream_vx / lift_dir_norm) + lift_yb = lift * (stream_vy / lift_dir_norm) + # Total lift force + R1, R2, R3 = lift_xb, lift_yb, 0 + # Total moment + M1, M2, M3 = -cpz * lift_yb, cpz * lift_xb, 0 + return R1, R2, R3, M1, M2, M3 diff --git a/rocketpy/rocket/aero_surface/fins/__init__.py b/rocketpy/rocket/aero_surface/fins/__init__.py index f1efc603a..941aa5465 100644 --- a/rocketpy/rocket/aero_surface/fins/__init__.py +++ b/rocketpy/rocket/aero_surface/fins/__init__.py @@ -1,3 +1,4 @@ from rocketpy.rocket.aero_surface.fins.elliptical_fins import EllipticalFins from rocketpy.rocket.aero_surface.fins.fins import Fins +from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins diff --git a/rocketpy/rocket/aero_surface/fins/fins.py b/rocketpy/rocket/aero_surface/fins/fins.py index 74f7922f7..e6cfb84ea 100644 --- a/rocketpy/rocket/aero_surface/fins/fins.py +++ b/rocketpy/rocket/aero_surface/fins/fins.py @@ -366,6 +366,66 @@ def fin_num_correction(n): else: return n / 2 + def compute_forces_and_moments( + self, + stream_velocity, + stream_speed, + stream_mach, + rho, + cp, + omega, + *args, + ): # pylint: disable=arguments-differ + """Computes the forces and moments acting on the aerodynamic surface. + + Parameters + ---------- + 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, + stream_mach, + rho, + cp, + ) + clf_delta, cld_omega, cant_angle_rad = self.roll_parameters + M3_forcing = ( + (1 / 2 * rho * stream_speed**2) + * self.reference_area + * self.reference_length + * clf_delta.get_value_opt(stream_mach) + * cant_angle_rad + ) + M3_damping = ( + (1 / 2 * rho * stream_speed) + * self.reference_area + * (self.reference_length) ** 2 + * cld_omega.get_value_opt(stream_mach) + * omega[2] + / 2 + ) + M3 = M3_forcing - M3_damping + return R1, R2, R3, M1, M2, M3 + 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. diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py new file mode 100644 index 000000000..1db2ef18d --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py @@ -0,0 +1,368 @@ +import warnings + +import numpy as np + +from rocketpy.plots.aero_surface_plots import _FreeFormFinsPlots +from rocketpy.prints.aero_surface_prints import _FreeFormFinsPrints + +from .fins import Fins + + +class FreeFormFins(Fins): + """Class that defines and holds information for a free form fin set. + + This class inherits from the Fins class. + + Note + ---- + Local coordinate system: + - Origin located at the top of the root chord. + - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + - Y axis perpendicular to the Z axis, in the span direction, positive upwards. + - X axis completes the right-handed coordinate system. + + See Also + -------- + Fins + + Attributes + ---------- + FreeFormFins.n : int + Number of fins in fin set. + FreeFormFins.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + FreeFormFins.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + FreeFormFins.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + FreeFormFins.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + FreeFormFins.root_chord : float + Fin root chord in meters. + FreeFormFins.span : float + Fin span in meters. + FreeFormFins.name : string + Name of fin set. + FreeFormFins.d : float + Reference diameter of the rocket, in meters. + FreeFormFins.ref_area : float + Reference area of the rocket, in m². + FreeFormFins.Af : float + Area of the longitudinal section of each fin in the set. + FreeFormFins.AR : float + Aspect ratio of each fin in the set + FreeFormFins.gamma_c : float + Fin mid-chord sweep angle. + FreeFormFins.Yma : float + Span wise position of the mean aerodynamic chord. + FreeFormFins.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + FreeFormFins.tau : float + Geometrical relation used to simplify lift and roll calculations. + FreeFormFins.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + FreeFormFins.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. + FreeFormFins.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + FreeFormFins.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + FreeFormFins.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + FreeFormFins.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. + FreeFormFins.clalpha : float + Lift coefficient slope. Has units of 1/rad. + FreeFormFins.mac_length : float + Mean aerodynamic chord length of the fin set. + FreeFormFins.mac_lead : float + Mean aerodynamic chord leading edge x coordinate. + """ + + def __init__( + self, + n, + shape_points, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Fins", + ): + """Initialize FreeFormFins class. + + Parameters + ---------- + n : int + Number of fins, from 2 to infinity. + shape_points : list + List of tuples (x, y) containing the coordinates of the fin's + geometry defining points. The point (0, 0) is the root leading edge. + Positive x is rearwards, positive y is upwards (span direction). + The shape will be interpolated between the points, in the order + they are given. The last point connects to the first point, and + represents the trailing edge. + 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. + 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 + """ + self.shape_points = shape_points + + down = False + for i in range(1, len(shape_points)): + if shape_points[i][1] > shape_points[i - 1][1] and down: + warnings.warn( + "Jagged fin shape detected. This may cause small inaccuracies " + "center of pressure and pitch moment calculations." + ) + break + if shape_points[i][1] < shape_points[i - 1][1]: + down = True + i += 1 + + root_chord = abs(shape_points[0][0] - shape_points[-1][0]) + ys = [point[1] for point in shape_points] + span = max(ys) - min(ys) + + 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 = _FreeFormFinsPrints(self) + self.plots = _FreeFormFinsPlots(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 = self.mac_lead + 0.25 * self.mac_length + 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 the fin set's geometrical parameters such as the + fin area, aspect ratio, and parameters related to roll movement. This + method uses the same calculations to those in OpenRocket for free-form + fin shapes. + + Returns + ------- + None + """ + # pylint: disable=invalid-name + # pylint: disable=too-many-locals + # Calculate the fin area (Af) using the Shoelace theorem (polygon area formula) + Af = 0 + for i in range(len(self.shape_points) - 1): + x1, y1 = self.shape_points[i] + x2, y2 = self.shape_points[i + 1] + Af += (y1 + y2) * (x1 - x2) + Af = abs(Af) / 2 + if Af < 1e-6: + raise ValueError("Fin area is too small. Check the shape_points.") + + # Calculate aspect ratio (AR) and lift interference factors + AR = 2 * self.span**2 / Af # Aspect ratio + tau = (self.span + self.rocket_radius) / self.rocket_radius + lift_interference_factor = 1 + 1 / tau + + # Calculate roll forcing interference factor using OpenRocket's approach + 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)) + ) + + # Define number of interpolation points along the span of the fin + points_per_line = 40 # Same as OpenRocket + + # Initialize arrays for leading/trailing edge and chord lengths + chord_lead = np.ones(points_per_line) * np.inf # Leading edge x coordinates + chord_trail = np.ones(points_per_line) * -np.inf # Trailing edge x coordinates + chord_length = np.zeros( + points_per_line + ) # Chord length for each spanwise section + + # Discretize fin shape and calculate chord length, leading, and trailing edges + for p in range(1, len(self.shape_points)): + x1, y1 = self.shape_points[p - 1] + x2, y2 = self.shape_points[p] + + # Compute corresponding points along the fin span (clamp to valid range) + prev_idx = int(y1 / self.span * (points_per_line - 1)) + curr_idx = int(y2 / self.span * (points_per_line - 1)) + prev_idx = np.clip(prev_idx, 0, points_per_line - 1) + curr_idx = np.clip(curr_idx, 0, points_per_line - 1) + + if prev_idx > curr_idx: + prev_idx, curr_idx = curr_idx, prev_idx + + # Compute intersection of fin edge with each spanwise section + for i in range(prev_idx, curr_idx + 1): + y = i * self.span / (points_per_line - 1) + if y1 != y2: + x = np.clip( + (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2, + min(x1, x2), + max(x1, x2), + ) + else: + x = x1 # Handle horizontal segments + + # Update leading and trailing edge positions + chord_lead[i] = min(chord_lead[i], x) + chord_trail[i] = max(chord_trail[i], x) + + # Update chord length + if y1 < y2: + chord_length[i] -= x + else: + chord_length[i] += x + + # Replace infinities and handle invalid values in chord_lead and chord_trail + for i in range(points_per_line): + if ( + np.isinf(chord_lead[i]) + or np.isinf(chord_trail[i]) + or np.isnan(chord_lead[i]) + or np.isnan(chord_trail[i]) + ): + chord_lead[i] = 0 + chord_trail[i] = 0 + if chord_length[i] < 0 or np.isnan(chord_length[i]): + chord_length[i] = 0 + if chord_length[i] > chord_trail[i] - chord_lead[i]: + chord_length[i] = chord_trail[i] - chord_lead[i] + + # Initialize integration variables for various aerodynamic and roll properties + radius = self.rocket_radius + total_area = 0 + mac_length = 0 # Mean aerodynamic chord length + mac_lead = 0 # Mean aerodynamic chord leading edge + mac_span = 0 # Mean aerodynamic chord spanwise position (Yma) + cos_gamma_sum = 0 # Sum of cosine of the sweep angle + roll_geometrical_constant = 0 + roll_damping_numerator = 0 + roll_damping_denominator = 0 + + # Perform integration over spanwise sections + dy = self.span / (points_per_line - 1) + for i in range(points_per_line): + chord = chord_trail[i] - chord_lead[i] + y = i * dy + + # Update integration variables + mac_length += chord * chord + mac_span += y * chord + mac_lead += chord_lead[i] * chord + total_area += chord + roll_geometrical_constant += chord_length[i] * (radius + y) ** 2 + roll_damping_numerator += radius**3 * chord / (radius + y) ** 2 + roll_damping_denominator += (radius + y) * chord + + # Update cosine of sweep angle (cos_gamma) + if i > 0: + dx = (chord_trail[i] + chord_lead[i]) / 2 - ( + chord_trail[i - 1] + chord_lead[i - 1] + ) / 2 + cos_gamma_sum += dy / np.hypot(dx, dy) + + # Finalize mean aerodynamic chord properties + mac_length *= dy + mac_span *= dy + mac_lead *= dy + total_area *= dy + roll_geometrical_constant *= dy + roll_damping_numerator *= dy + roll_damping_denominator *= dy + + mac_length /= total_area + mac_span /= total_area + mac_lead /= total_area + cos_gamma = cos_gamma_sum / (points_per_line - 1) + + # Store computed values + self.Af = Af # Fin area + self.AR = AR # Aspect ratio + self.gamma_c = np.arccos(cos_gamma) # Sweep angle + self.Yma = mac_span # Mean aerodynamic chord spanwise position + self.mac_length = mac_length + self.mac_lead = mac_lead + self.tau = tau + self.roll_geometrical_constant = roll_geometrical_constant + self.lift_interference_factor = lift_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + self.roll_damping_interference_factor = 1 + ( + roll_damping_numerator / roll_damping_denominator + ) + + # Evaluate the shape and finalize geometry + self.evaluate_shape() + + def evaluate_shape(self): + x_array, y_array = zip(*self.shape_points) + self.shape_vec = [np.array(x_array), np.array(y_array)] + + def info(self): + self.prints.geometry() + self.prints.lift() + + def all_info(self): + self.prints.all() + self.plots.all() diff --git a/rocketpy/rocket/aero_surface/generic_surface.py b/rocketpy/rocket/aero_surface/generic_surface.py new file mode 100644 index 000000000..1a0e1015d --- /dev/null +++ b/rocketpy/rocket/aero_surface/generic_surface.py @@ -0,0 +1,445 @@ +import copy +import csv +import math + +import numpy as np + +from rocketpy.mathutils import Function +from rocketpy.mathutils.vector_matrix import Matrix, Vector + + +class GenericSurface: + """Defines a generic aerodynamic surface with custom force and moment + coefficients. The coefficients can be nonlinear functions of the angle of + attack, sideslip angle, Mach number, Reynolds number, pitch rate, yaw rate + and roll rate.""" + + def __init__( + self, + reference_area, + reference_length, + coefficients, + center_of_pressure=(0, 0, 0), + name="Generic Surface", + ): + """Create a generic aerodynamic surface, defined by its aerodynamic + coefficients. This surface is used to model any aerodynamic surface + that does not fit the predefined classes. + + Important + --------- + All the aerodynamic coefficients can be input as callable functions of + angle of attack, angle of sideslip, Mach number, Reynolds number, + pitch rate, yaw rate and roll rate. For CSV files, the header must + contain at least one of the following: "alpha", "beta", "mach", + "reynolds", "pitch_rate", "yaw_rate" and "roll_rate". + + See Also + -------- + :ref:`genericsurfaces`. + + Parameters + ---------- + reference_area : int, float + Reference area of the aerodynamic surface. Has the unit of meters + squared. Commonly defined as the rocket's cross-sectional area. + reference_length : int, float + Reference length of the aerodynamic surface. Has the unit of meters. + Commonly defined as the rocket's diameter. + 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 + aerodynamic surface. The default value is (0, 0, 0). + 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 + 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 + + 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, + rho, + stream_speed, + alpha, + beta, + mach, + reynolds, + pitch_rate, + yaw_rate, + roll_rate, + ): + """Compute the aerodynamic forces and moments from the aerodynamic + coefficients. + + Parameters + ---------- + rho : float + Air density. + stream_speed : float + Magnitude of the airflow speed. + alpha : float + Angle of attack in radians. + beta : float + Sideslip angle in radians. + mach : float + Mach number. + reynolds : float + Reynolds number. + pitch_rate : float + Pitch rate in radians per second. + yaw_rate : float + Yaw rate in radians per second. + roll_rate : float + Roll rate in radians per second. + + Returns + ------- + tuple of float + The aerodynamic forces (lift, side_force, drag) and moments + (pitch, yaw, roll) in the body frame. + """ + # Precompute common values + dyn_pressure_area = 0.5 * rho * stream_speed**2 * self.reference_area + dyn_pressure_area_length = dyn_pressure_area * self.reference_length + + # Compute aerodynamic forces + lift = dyn_pressure_area * self.cL( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + side = dyn_pressure_area * self.cQ( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + drag = dyn_pressure_area * self.cD( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + # Compute aerodynamic moments + pitch = dyn_pressure_area_length * self.cm( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + yaw = dyn_pressure_area_length * self.cn( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + roll = dyn_pressure_area_length * self.cl( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + return lift, side, drag, pitch, yaw, roll + + def compute_forces_and_moments( + self, + stream_velocity, + stream_speed, + stream_mach, + rho, + cp, + omega, + reynolds, + ): + """Computes the forces and moments acting on the aerodynamic surface. + Used in each time step of the simulation. This method is valid for + both linear and nonlinear aerodynamic coefficients. + + Parameters + ---------- + 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. + reynolds : float + Reynolds number. + omega: tuple of 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. + """ + # Stream velocity in standard aerodynamic frame + stream_velocity = -stream_velocity + + # Angles of attack and sideslip + alpha = np.arctan2(stream_velocity[1], stream_velocity[2]) + beta = np.arctan2(-stream_velocity[0], stream_velocity[2]) + + # Compute aerodynamic forces and moments + lift, side, drag, pitch, yaw, roll = self._compute_from_coefficients( + rho, + stream_speed, + alpha, + beta, + stream_mach, + reynolds, + omega[0], + omega[1], + omega[2], + ) + + # Conversion from aerodynamic frame to body frame + rotation_matrix = Matrix( + [ + [1, 0, 0], + [0, math.cos(alpha), -math.sin(alpha)], + [0, math.sin(alpha), math.cos(alpha)], + ] + ) @ Matrix( + [ + [math.cos(beta), 0, -math.sin(beta)], + [0, 1, 0], + [math.sin(beta), 0, math.cos(beta)], + ] + ) + R1, R2, R3 = rotation_matrix @ Vector([side, -lift, -drag]) + + # Dislocation of the aerodynamic application point to CDM + M1, M2, M3 = Vector([pitch, yaw, roll]) + (cp ^ Vector([R1, R2, R3])) + + return R1, R2, R3, M1, M2, M3 + + def _process_input(self, input_data, coeff_name): + """Process the input data, either as a CSV file or a callable function. + + Parameters + ---------- + input_data : str or callable + Input data to be processed, either a path to a CSV or a callable. + coeff_name : str + Name of the coefficient being processed for error reporting. + + Returns + ------- + Function + Function object with 7 input arguments (alpha, beta, mach, reynolds, + pitch_rate, yaw_rate, roll_rate). + """ + if isinstance(input_data, str): + # Input is assumed to be a file path to a CSV + return self.__load_csv(input_data, coeff_name) + elif isinstance(input_data, Function): + if input_data.__dom_dim__ != 7: + raise ValueError( + f"{coeff_name} function must have 7 input arguments" + " (alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate)." + ) + return input_data + elif callable(input_data): + # Check if callable has 7 inputs (alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + if input_data.__code__.co_argcount != 7: + raise ValueError( + f"{coeff_name} function must have 7 input arguments" + " (alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate)." + ) + return input_data + elif input_data == 0: + return Function( + lambda alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate: 0, + [ + 'alpha', + 'beta', + 'mach', + 'reynolds', + 'pitch_rate', + 'yaw_rate', + 'roll_rate', + ], + [coeff_name], + interpolation='linear', + extrapolation='natural', + ) + else: + raise TypeError( + f"Invalid input for {coeff_name}: must be a CSV file path" + " or a callable." + ) + + def __load_csv(self, file_path, coeff_name): + """Load a CSV file and create a Function object with the correct number + of arguments. The CSV file must have a header that specifies the + independent variables that are used. + + Parameters + ---------- + file_path : str + Path to the CSV file. + coeff_name : str + Name of the coefficient being processed. + + Returns + ------- + Function + Function object with 7 input arguments (alpha, beta, mach, reynolds, + pitch_rate, yaw_rate, roll_rate). + """ + try: + with open(file_path, mode='r') as file: + reader = csv.reader(file) + header = next(reader) + except (FileNotFoundError, IOError) as 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}.") + + # TODO make header strings flexible (e.g. 'alpha', 'Alpha', 'ALPHA') + independent_vars = [ + 'alpha', + 'beta', + 'mach', + 'reynolds', + 'pitch_rate', + 'yaw_rate', + 'roll_rate', + ] + present_columns = [col for col in independent_vars if col in header] + + # Check that the last column is not an independent variable + if header[-1] in independent_vars: + raise ValueError( + f"Last column in {coeff_name} CSV must be the coefficient" + " value, not an independent variable." + ) + + # Ensure that at least one independent variable is present + if not present_columns: + raise ValueError(f"No independent variables found in {coeff_name} CSV.") + + # Initialize the CSV-based function + csv_func = Function( + file_path, + interpolation='linear', + extrapolation='natural', + ) + + # Create a mask for the presence of each independent variable + # save on self to avoid loss of scope + _mask = [1 if col in present_columns else 0 for col in independent_vars] + + # Generate a lambda that applies only the relevant arguments to csv_func + def wrapper(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate): + args = [alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate] + # Select arguments that correspond to present variables + selected_args = [arg for arg, m in zip(args, _mask) if m] + return csv_func(*selected_args) + + # Create the interpolation function + func = Function( + wrapper, + independent_vars, + [coeff_name], + interpolation='linear', + extrapolation='natural', + ) + return func diff --git a/rocketpy/rocket/aero_surface/linear_generic_surface.py b/rocketpy/rocket/aero_surface/linear_generic_surface.py new file mode 100644 index 000000000..fba2b1db8 --- /dev/null +++ b/rocketpy/rocket/aero_surface/linear_generic_surface.py @@ -0,0 +1,410 @@ +from rocketpy.mathutils import Function +from rocketpy.plots.aero_surface_plots import _LinearGenericSurfacePlots +from rocketpy.prints.aero_surface_prints import _LinearGenericSurfacePrints +from rocketpy.rocket.aero_surface.generic_surface import GenericSurface + + +class LinearGenericSurface(GenericSurface): + """Class that defines a generic linear aerodynamic surface. This class is + used to define aerodynamic surfaces that have aerodynamic coefficients + defined as linear functions of the coefficients derivatives.""" + + def __init__( + self, + reference_area, + reference_length, + coefficients, + center_of_pressure=(0, 0, 0), + name="Generic Linear Surface", + ): + """Create a generic linear aerodynamic surface, defined by its + aerodynamic coefficients derivatives. This surface is used to model any + aerodynamic surface that does not fit the predefined classes. + + Important + --------- + All the aerodynamic coefficients can be input as callable functions of + angle of attack, angle of sideslip, Mach number, Reynolds number, + pitch rate, yaw rate and roll rate. For CSV files, the header must + contain at least one of the following: "alpha", "beta", "mach", + "reynolds", "pitch_rate", "yaw_rate" and "roll_rate". + + See Also + -------- + :ref:`genericsurfaces`. + + Parameters + ---------- + reference_area : int, float + Reference area of the aerodynamic surface. Has the unit of meters + squared. Commonly defined as the rocket's cross-sectional area. + reference_length : int, float + Reference length of the aerodynamic surface. Has the unit of meters. + Commonly defined as the rocket's diameter. + 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 + aerodynamic surface. The default value is (0, 0, 0). + name : str + Name of the aerodynamic surface. Default is 'GenericSurface'. + """ + + super().__init__( + reference_area=reference_area, + reference_length=reference_length, + coefficients=coefficients, + center_of_pressure=center_of_pressure, + name=name, + ) + + self.compute_all_coefficients() + + self.prints = _LinearGenericSurfacePrints(self) + self.plots = _LinearGenericSurfacePlots(self) + + 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.""" + + def total_coefficient( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ): + return ( + c_0(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + + c_alpha(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + * alpha + + c_beta(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + * beta + ) + + return Function( + total_coefficient, + [ + 'alpha', + 'beta', + 'mach', + 'reynolds', + 'pitch_rate', + 'yaw_rate', + 'roll_rate', + ], + ['coefficient'], + ) + + def compute_damping_coefficient(self, c_p, c_q, c_r): + """Compute the damping coefficient from the derivatives of the + aerodynamic coefficients.""" + + def total_coefficient( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ): + return ( + c_p(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + * pitch_rate + + c_q(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + * yaw_rate + + c_r(alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate) + * roll_rate + ) + + return Function( + total_coefficient, + [ + 'alpha', + 'beta', + 'mach', + 'reynolds', + 'pitch_rate', + 'yaw_rate', + 'roll_rate', + ], + ['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 + ) + self.cLd = self.compute_damping_coefficient(self.cL_p, self.cL_q, self.cL_r) + + self.cQf = self.compute_forcing_coefficient( + self.cQ_0, self.cQ_alpha, self.cQ_beta + ) + self.cQd = self.compute_damping_coefficient(self.cQ_p, self.cQ_q, self.cQ_r) + + self.cDf = self.compute_forcing_coefficient( + self.cD_0, self.cD_alpha, self.cD_beta + ) + self.cDd = self.compute_damping_coefficient(self.cD_p, self.cD_q, self.cD_r) + + self.cmf = self.compute_forcing_coefficient( + self.cm_0, self.cm_alpha, self.cm_beta + ) + self.cmd = self.compute_damping_coefficient(self.cm_p, self.cm_q, self.cm_r) + + self.cnf = self.compute_forcing_coefficient( + self.cn_0, self.cn_alpha, self.cn_beta + ) + self.cnd = self.compute_damping_coefficient(self.cn_p, self.cn_q, self.cn_r) + + self.clf = self.compute_forcing_coefficient( + self.cl_0, self.cl_alpha, self.cl_beta + ) + self.cld = self.compute_damping_coefficient(self.cl_p, self.cl_q, self.cl_r) + + def _compute_from_coefficients( + self, + rho, + stream_speed, + alpha, + beta, + mach, + reynolds, + pitch_rate, + yaw_rate, + roll_rate, + ): + """Compute the aerodynamic forces and moments from the aerodynamic + coefficients. + + Parameters + ---------- + rho : float + Air density. + stream_speed : float + Magnitude of the airflow speed. + alpha : float + Angle of attack in radians. + beta : float + Sideslip angle in radians. + mach : float + Mach number. + reynolds : float + Reynolds number. + pitch_rate : float + Pitch rate in radians per second. + yaw_rate : float + Yaw rate in radians per second. + roll_rate : float + Roll rate in radians per second. + + Returns + ------- + tuple of float + The aerodynamic forces (lift, side_force, drag) and moments + (pitch, yaw, roll) in the body frame. + """ + # Precompute common values + dyn_pressure_area = 0.5 * rho * stream_speed**2 * self.reference_area + dyn_pressure_area_damping = ( + 0.5 * rho * stream_speed * self.reference_area * self.reference_length / 2 + ) + dyn_pressure_area_length = dyn_pressure_area * self.reference_length + dyn_pressure_area_length_damping = ( + 0.5 + * rho + * stream_speed + * self.reference_area + * self.reference_length**2 + / 2 + ) + + # Compute aerodynamic forces + lift = dyn_pressure_area * self.cLf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_damping * self.cLd( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + side = dyn_pressure_area * self.cQf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_damping * self.cQd( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + drag = dyn_pressure_area * self.cDf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_damping * self.cDd( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + # Compute aerodynamic moments + pitch = dyn_pressure_area_length * self.cmf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_length_damping * self.cmd( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + yaw = dyn_pressure_area_length * self.cnf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_length_damping * self.cnd( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + roll = dyn_pressure_area_length * self.clf( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) - dyn_pressure_area_length_damping * self.cld( + alpha, beta, mach, reynolds, pitch_rate, yaw_rate, roll_rate + ) + + return lift, side, drag, pitch, yaw, roll diff --git a/rocketpy/rocket/aero_surface/nose_cone.py b/rocketpy/rocket/aero_surface/nose_cone.py index 8d98113ae..7d59473e3 100644 --- a/rocketpy/rocket/aero_surface/nose_cone.py +++ b/rocketpy/rocket/aero_surface/nose_cone.py @@ -446,10 +446,6 @@ def evaluate_lift_coefficient(self): 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.radius_ratio**2, "Mach", diff --git a/rocketpy/rocket/aero_surface/tail.py b/rocketpy/rocket/aero_surface/tail.py index 5785232fd..4f2783e02 100644 --- a/rocketpy/rocket/aero_surface/tail.py +++ b/rocketpy/rocket/aero_surface/tail.py @@ -163,10 +163,6 @@ def evaluate_lift_coefficient(self): 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 * ( diff --git a/rocketpy/rocket/components.py b/rocketpy/rocket/components.py index 2d580de7e..91bd05fa0 100644 --- a/rocketpy/rocket/components.py +++ b/rocketpy/rocket/components.py @@ -180,7 +180,7 @@ def clear(self): self._components.clear() def sort_by_position(self, reverse=False): - """Sort the list of components by position. + """Sort the list of components by z axis position. Parameters ---------- @@ -192,4 +192,4 @@ def sort_by_position(self, reverse=False): ------- None """ - self._components.sort(key=lambda x: x.position, reverse=reverse) + self._components.sort(key=lambda x: x.position.z, reverse=reverse) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index a5964808f..ec3bfc0ec 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1,3 +1,4 @@ +import math import warnings import numpy as np @@ -17,6 +18,8 @@ Tail, TrapezoidalFins, ) +from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins +from rocketpy.rocket.aero_surface.generic_surface import GenericSurface from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute from rocketpy.tools import parallel_axis_theorem_from_com @@ -107,6 +110,11 @@ class Rocket: Rocket.aerodynamic_surfaces : list Collection of aerodynamic surfaces of the rocket. Holds Nose cones, Fin sets, and Tails. + Rocket.surfaces_cp_to_cdm : dict + Dictionary containing the relative position of each aerodynamic surface + center of pressure to the rocket's center of mass. The key is the + aerodynamic surface object and the value is the relative position Vector + in meters. Rocket.parachutes : list Collection of parachutes of the rocket. Rocket.air_brakes : list @@ -290,6 +298,8 @@ def __init__( # pylint: disable=too-many-statements self.area = np.pi * self.radius**2 # Eccentricity data initialization + self.cm_eccentricity_x = 0 + self.cm_eccentricity_y = 0 self.cp_eccentricity_x = 0 self.cp_eccentricity_y = 0 self.thrust_eccentricity_y = 0 @@ -301,6 +311,7 @@ def __init__( # pylint: disable=too-many-statements self.air_brakes = [] self.sensors = Components() self.aerodynamic_surfaces = Components() + self.surfaces_cp_to_cdm = {} self.rail_buttons = Components() self.cp_position = Function( @@ -525,18 +536,50 @@ def evaluate_center_of_pressure(self): # Calculate total lift coefficient derivative and center of pressure if len(self.aerodynamic_surfaces) > 0: for aero_surface, position in self.aerodynamic_surfaces: + if isinstance(aero_surface, GenericSurface): + continue # ref_factor corrects lift for different reference areas ref_factor = (aero_surface.rocket_radius / self.radius) ** 2 self.total_lift_coeff_der += ref_factor * aero_surface.clalpha self.cp_position += ( ref_factor * aero_surface.clalpha - * (position - self._csys * aero_surface.cpz) + * (position.z - self._csys * aero_surface.cpz) ) - self.cp_position /= self.total_lift_coeff_der - + # Avoid errors when only generic surfaces are added + if self.total_lift_coeff_der.get_value(0) != 0: + self.cp_position /= self.total_lift_coeff_der return self.cp_position + def evaluate_surfaces_cp_to_cdm(self): + """Calculates the relative position of each aerodynamic surface center + of pressure to the rocket's center of dry mass in Body Axes Coordinate + System. + + Returns + ------- + self.surfaces_cp_to_cdm : dict + Dictionary mapping the relative position of each aerodynamic + surface center of pressure to the rocket's center of mass. + """ + for surface, position in self.aerodynamic_surfaces: + self.evaluate_single_surface_cp_to_cdm(surface, position) + return self.surfaces_cp_to_cdm + + def evaluate_single_surface_cp_to_cdm(self, surface, position): + """Calculates the relative position of each aerodynamic surface + center of pressure to the rocket's center of dry mass in Body Axes + Coordinate System.""" + pos = Vector( + [ + (position.x - self.cm_eccentricity_x) * self._csys - surface.cpx, + (position.y - self.cm_eccentricity_y) - surface.cpy, + (position.z - self.center_of_dry_mass_position) * self._csys + - surface.cpz, + ] + ) + self.surfaces_cp_to_cdm[surface] = pos + def evaluate_stability_margin(self): """Calculates the stability margin of the rocket as a function of mach number and time. @@ -917,6 +960,7 @@ def add_motor(self, motor, position): # pylint: disable=too-many-statements self.evaluate_reduced_mass() self.evaluate_thrust_to_weight() self.evaluate_center_of_pressure() + self.evaluate_surfaces_cp_to_cdm() self.evaluate_stability_margin() self.evaluate_static_margin() self.evaluate_com_to_cdm_function() @@ -932,7 +976,7 @@ def add_surfaces(self, surfaces, positions): surfaces : list, AeroSurface, NoseCone, TrapezoidalFins, EllipticalFins, Tail Aerodynamic surface to be added to the rocket. Can be a list of AeroSurface if more than one surface is to be added. - positions : int, float, list + positions : int, float, list, tuple, Vector Position, in m, of the aerodynamic surface's center of pressure relative to the user defined rocket coordinate system. If a list is passed, it will correspond to the position of each item @@ -952,16 +996,21 @@ def add_surfaces(self, surfaces, positions): ------- None """ - if not isinstance(surfaces, list): - surfaces = [surfaces] - positions = [positions] - - for surface, position in zip(surfaces, positions): - if isinstance(surface, RailButtons): - surface.rocket_radius = surface.rocket_radius or self.radius - self.rail_buttons.add(surface, position) - else: + try: + for surface, position in zip(surfaces, positions): + if not isinstance(position, (Vector, tuple, list)): + position = Vector([0, 0, position]) + else: + position = Vector(position) self.aerodynamic_surfaces.add(surface, position) + self.evaluate_single_surface_cp_to_cdm(surface, position) + except TypeError: + if not isinstance(positions, (Vector, tuple, list)): + positions = Vector([0, 0, positions]) + else: + positions = Vector(positions) + self.aerodynamic_surfaces.add(surfaces, positions) + self.evaluate_single_surface_cp_to_cdm(surfaces, positions) self.evaluate_center_of_pressure() self.evaluate_stability_margin() @@ -1274,6 +1323,84 @@ def add_elliptical_fins( self.add_surfaces(fin_set, position) return fin_set + def add_free_form_fins( + self, + n, + shape_points, + position, + cant_angle=0.0, + radius=None, + airfoil=None, + name="Fins", + ): + """Create a free form fin set, storing its parameters as part of the + aerodynamic_surfaces list. Its parameters are the axial position along + the rocket and its derivative of the coefficient of lift in respect to + angle of attack. + + Parameters + ---------- + n : int + Number of fins, from 2 to infinity. + shape_points : list + List of tuples (x, y) containing the coordinates of the fin's + geometry defining points. The point (0, 0) is the root leading edge. + Positive x is rearwards, positive y is upwards (span direction). + The shape will be interpolated between the points, in the order + they are given. The last point connects to the first point. + position : int, float + Fin set position relative to the rocket's coordinate system. + By fin set position, understand the point belonging to the root + chord which is highest in the rocket coordinate system (i.e. + the point closest to the nose cone tip). + + See Also + -------- + :ref:`positions` + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + radius : int, float, optional + Reference fuselage radius where the fins are located. This is used + to calculate lift coefficient and to draw the rocket. If None, + which is default, the rocket radius will be used. + 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". + + Returns + ------- + fin_set : FreeFormFins + Fin set object created. + """ + + # Modify radius if not given, use rocket radius, otherwise use given. + radius = radius if radius is not None else self.radius + + fin_set = FreeFormFins( + n, + shape_points, + radius, + cant_angle, + airfoil, + name, + ) + + # Add fin set to the list of aerodynamic surfaces + self.add_surfaces(fin_set, position) + return fin_set + def add_parachute( self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0) ): @@ -1564,13 +1691,22 @@ def set_rail_buttons( rail_buttons : RailButtons RailButtons object created """ + radius = radius or self.radius buttons_distance = abs(upper_button_position - lower_button_position) rail_buttons = RailButtons( - buttons_distance=buttons_distance, angular_position=angular_position + buttons_distance=buttons_distance, + angular_position=angular_position, + rocket_radius=radius, ) self.rail_buttons = Components() - rail_buttons.rocket_radius = radius or self.radius - self.rail_buttons.add(rail_buttons, lower_button_position) + position = Vector( + [ + radius * -math.sin(math.radians(angular_position)), + radius * math.cos(math.radians(angular_position)), + lower_button_position, + ] + ) + self.rail_buttons.add(rail_buttons, position) return rail_buttons def add_cm_eccentricity(self, x, y): @@ -1583,25 +1719,31 @@ def add_cm_eccentricity(self, x, y): ---------- x : float Distance in meters by which the CM is to be translated in - the x direction relative to geometrical center line. + the x direction relative to geometrical center line. The x axis + is defined according to the body axes coordinate system. y : float Distance in meters by which the CM is to be translated in - the y direction relative to geometrical center line. + the y direction relative to geometrical center line. The y axis + is defined according to the body axes coordinate system. Returns ------- self : Rocket Object of the Rocket class. + See Also + -------- + :ref:`rocketaxes` + Notes ----- Should not be used together with add_cp_eccentricity and add_thrust_eccentricity. """ - self.cp_eccentricity_x = -x - self.cp_eccentricity_y = -y - self.thrust_eccentricity_y = -x - self.thrust_eccentricity_x = -y + self.cm_eccentricity_x = x + self.cm_eccentricity_y = y + self.add_cp_eccentricity(-x, -y) + self.add_thrust_eccentricity(-x, -y) return self def add_cp_eccentricity(self, x, y): @@ -1614,14 +1756,22 @@ def add_cp_eccentricity(self, x, y): x : float Distance in meters by which the CP is to be translated in the x direction relative to the center of mass axial line. + The x axis is defined according to the body axes coordinate + system. y : float Distance in meters by which the CP is to be translated in the y direction relative to the center of mass axial line. + The y axis is defined according to the body axes coordinate + system. Returns ------- self : Rocket Object of the Rocket class. + + See Also + -------- + :ref:`rocketaxes` """ self.cp_eccentricity_x = x self.cp_eccentricity_y = y @@ -1636,16 +1786,22 @@ def add_thrust_eccentricity(self, x, y): x : float Distance in meters by which the line of action of the thrust force is to be translated in the x direction - relative to the center of mass axial line. + relative to the center of mass axial line. The x axis + is defined according to the body axes coordinate system. y : float Distance in meters by which the line of action of the thrust force is to be translated in the x direction - relative to the center of mass axial line. + relative to the center of mass axial line. The y axis + is defined according to the body axes coordinate system. Returns ------- self : Rocket Object of the Rocket class. + + See Also + -------- + :ref:`rocketaxes` """ self.thrust_eccentricity_y = x self.thrust_eccentricity_x = y diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index e7b53754a..70ddff6c1 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -736,10 +736,11 @@ def __simulate(self, verbose): for parachute in node.parachutes: # Calculate and save pressure signal - noisy_pressure, height_above_ground_level = ( - self.__calculate_and_save_pressure_signals( - parachute, node.t, self.y_sol[2] - ) + ( + noisy_pressure, + height_above_ground_level, + ) = self.__calculate_and_save_pressure_signals( + parachute, node.t, self.y_sol[2] ) if parachute.triggerfunc( noisy_pressure, @@ -975,12 +976,13 @@ def __simulate(self, verbose): ) for parachute in overshootable_node.parachutes: # Calculate and save pressure signal - noisy_pressure, height_above_ground_level = ( - self.__calculate_and_save_pressure_signals( - parachute, - overshootable_node.t, - overshootable_node.y_sol[2], - ) + ( + noisy_pressure, + height_above_ground_level, + ) = self.__calculate_and_save_pressure_signals( + parachute, + overshootable_node.t, + overshootable_node.y_sol[2], ) # Check for parachute trigger @@ -1230,7 +1232,7 @@ def effective_1rl(self): rail_buttons = self.rocket.rail_buttons[0] upper_r_button = ( rail_buttons.component.buttons_distance * self.rocket._csys - + rail_buttons.position + + rail_buttons.position.z ) except IndexError: # No rail buttons defined upper_r_button = nozzle @@ -1245,7 +1247,7 @@ def effective_2rl(self): nozzle = self.rocket.nozzle_position try: rail_buttons = self.rocket.rail_buttons[0] - lower_r_button = rail_buttons.position + lower_r_button = rail_buttons.position.z except IndexError: # No rail buttons defined lower_r_button = nozzle effective_2rl = self.rail_length - abs(nozzle - lower_r_button) @@ -1464,16 +1466,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 +1511,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 +1638,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 @@ -1788,13 +1763,9 @@ def u_dot_generalized( # Get rocket velocity in body frame velocity_in_body_frame = Kt @ v # Calculate lift and moment for each component of the rocket - for aero_surface, position in self.rocket.aerodynamic_surfaces: - comp_cpz = ( - position - self.rocket.center_of_dry_mass_position - ) * self.rocket._csys - aero_surface.cpz - comp_cp = Vector([0, 0, comp_cpz]) - reference_area = aero_surface.reference_area - reference_length = aero_surface.reference_length + 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_vb = velocity_in_body_frame + (w ^ comp_cp) # Wind velocity at component altitude @@ -1804,55 +1775,33 @@ def u_dot_generalized( # Component freestream velocity in body frame comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) comp_stream_velocity = comp_wind_vb - comp_vb - comp_stream_vx_b, comp_stream_vy_b, comp_stream_vz_b = comp_stream_velocity comp_stream_speed = abs(comp_stream_velocity) comp_stream_mach = comp_stream_speed / speed_of_sound - # 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: # TODO: maybe try/except - # 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, comp_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_cpz) * comp_lift_yb - M2 += (comp_cpz) * comp_lift_xb - # Calculates Roll Moment - try: - clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters - M3_forcing = ( - (1 / 2 * rho * comp_stream_speed**2) - * reference_area - * reference_length - * clf_delta.get_value_opt(comp_stream_mach) - * cant_angle_rad - ) - M3_damping = ( - (1 / 2 * rho * comp_stream_speed) - * reference_area - * (reference_length) ** 2 - * cld_omega.get_value_opt(comp_stream_mach) - * omega3 - / 2 - ) - M3 += M3_forcing - M3_damping - except AttributeError: - pass + # 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 thrust = self.rocket.motor.thrust.get_value_opt(t) @@ -2058,37 +2007,44 @@ def get_solution_at_time(self, t, atol=1e-3): # Transform solution array into Functions @funcify_method("Time (s)", "X (m)", "spline", "constant") def x(self): - """Rocket x position as a Function of time.""" + """Rocket x position relative to the launch pad as a Function of + time.""" return self.solution_array[:, [0, 1]] @funcify_method("Time (s)", "Y (m)", "spline", "constant") def y(self): - """Rocket y position as a Function of time.""" + """Rocket y position relative to the lauch pad as a Function of + time.""" return self.solution_array[:, [0, 2]] @funcify_method("Time (s)", "Z (m)", "spline", "constant") def z(self): - """Rocket z position as a Function of time.""" + """Rocket z position relative to the launch pad as a Function of + time.""" return self.solution_array[:, [0, 3]] @funcify_method("Time (s)", "Altitude AGL (m)", "spline", "constant") def altitude(self): - """Rocket altitude above ground level as a Function of time.""" + """Rocket altitude above ground level as a Function of time. Ground + level is defined by the environment elevation.""" return self.z - self.env.elevation @funcify_method("Time (s)", "Vx (m/s)", "spline", "zero") def vx(self): - """Rocket x velocity as a Function of time.""" + """Velocity of the rocket center of dry mass in the direction of + the rocket x-axis as a Function of time.""" return self.solution_array[:, [0, 4]] @funcify_method("Time (s)", "Vy (m/s)", "spline", "zero") def vy(self): - """Rocket y velocity as a Function of time.""" + """Velocity of the rocket center of dry mass in the direction of + the rocket y-axis as a Function of time.""" return self.solution_array[:, [0, 5]] @funcify_method("Time (s)", "Vz (m/s)", "spline", "zero") def vz(self): - """Rocket z velocity as a Function of time.""" + """Velocity of the rocket center of dry mass in the direction of + the rocket z-axis as a Function of time.""" return self.solution_array[:, [0, 6]] @funcify_method("Time (s)", "e0", "spline", "constant") @@ -2113,88 +2069,94 @@ def e3(self): @funcify_method("Time (s)", "ω1 (rad/s)", "spline", "zero") def w1(self): - """Rocket angular velocity ω1 as a Function of time.""" + """Angular velocity of the rocket in the x-axis as a Function of time. + Sometimes referred to as pitch rate (q).""" return self.solution_array[:, [0, 11]] @funcify_method("Time (s)", "ω2 (rad/s)", "spline", "zero") def w2(self): - """Rocket angular velocity ω2 as a Function of time.""" + """Angular velocity of the rocket in the y-axis as a Function of time. + Sometimes referred to as yaw rate (r).""" return self.solution_array[:, [0, 12]] @funcify_method("Time (s)", "ω3 (rad/s)", "spline", "zero") def w3(self): - """Rocket angular velocity ω3 as a Function of time.""" + """Angular velocity of the rocket in the z-axis as a Function of time. + Sometimes referred to as roll rate (p).""" return self.solution_array[:, [0, 13]] # Process second type of outputs - accelerations components @funcify_method("Time (s)", "Ax (m/s²)", "spline", "zero") def ax(self): - """Rocket x acceleration as a Function of time.""" + """Acceleration of the rocket center of dry mass in the direction of + the rocket x-axis as a Function of time.""" return self.__evaluate_post_process[:, [0, 1]] @funcify_method("Time (s)", "Ay (m/s²)", "spline", "zero") def ay(self): - """Rocket y acceleration as a Function of time.""" + """Acceleration of the rocket center of dry mass in the direction of + the rocket y-axis as a Function of time.""" return self.__evaluate_post_process[:, [0, 2]] @funcify_method("Time (s)", "Az (m/s²)", "spline", "zero") def az(self): - """Rocket z acceleration as a Function of time.""" + """Acceleration of the rocket center of dry mass in the direction of + the rocket z-axis as a Function of time.""" return self.__evaluate_post_process[:, [0, 3]] @funcify_method("Time (s)", "α1 (rad/s²)", "spline", "zero") def alpha1(self): - """Rocket angular acceleration α1 as a Function of time.""" + """Angular acceleration of the rocket in the x-axis as a Function of + time. Sometimes referred to as pitch acceleration.""" return self.__evaluate_post_process[:, [0, 4]] @funcify_method("Time (s)", "α2 (rad/s²)", "spline", "zero") def alpha2(self): - """Rocket angular acceleration α2 as a Function of time.""" + """Angular acceleration of the rocket in the y-axis as a Function of + time. Sometimes referred to as yaw acceleration.""" return self.__evaluate_post_process[:, [0, 5]] @funcify_method("Time (s)", "α3 (rad/s²)", "spline", "zero") def alpha3(self): - """Rocket angular acceleration α3 as a Function of time.""" + """Angular acceleration of the rocket in the z-axis as a Function of + time. Sometimes referred to as roll acceleration.""" return self.__evaluate_post_process[:, [0, 6]] # Process third type of outputs - Temporary values @funcify_method("Time (s)", "R1 (N)", "spline", "zero") def R1(self): - """Aerodynamic force along the first axis that is perpendicular to the - rocket's axis of symmetry as a Function of time.""" + """Aerodynamic force along the x-axis of the rocket as a Function of + time.""" return self.__evaluate_post_process[:, [0, 7]] @funcify_method("Time (s)", "R2 (N)", "spline", "zero") def R2(self): - """Aerodynamic force along the second axis that is perpendicular to the - rocket's axis of symmetry as a Function of time.""" + """Aerodynamic force along the y-axis of the rocket as a Function of + time.""" return self.__evaluate_post_process[:, [0, 8]] @funcify_method("Time (s)", "R3 (N)", "spline", "zero") def R3(self): - """Aerodynamic force along the rocket's axis of symmetry as a + """Aerodynamic force along the z-axis (rocket's axis of symmetry) as a Function of time.""" return self.__evaluate_post_process[:, [0, 9]] - @funcify_method("Time (s)", "M1 (Nm)", "spline", "zero") + @funcify_method("Time (s)", "M1 (Nm)", "linear", "zero") def M1(self): - """Aerodynamic bending moment in the same direction as the axis that is - perpendicular to the rocket's axis of symmetry as a Function of - time. - """ + """Aerodynamic moment in the rocket x-axis as a Function of time. + Sometimes referred to as pitch moment.""" return self.__evaluate_post_process[:, [0, 10]] - @funcify_method("Time (s)", "M2 (Nm)", "spline", "zero") + @funcify_method("Time (s)", "M2 (Nm)", "linear", "zero") def M2(self): - """Aerodynamic bending moment in the same direction as the axis that is - perpendicular to the rocket's axis of symmetry as a Function - of time.""" + """Aerodynamic moment in the rocket y-axis as a Function of time. + Sometimes referred to as yaw moment.""" return self.__evaluate_post_process[:, [0, 11]] - @funcify_method("Time (s)", "M3 (Nm)", "spline", "zero") + @funcify_method("Time (s)", "M3 (Nm)", "linear", "zero") def M3(self): - """Aerodynamic bending moment in the same direction as the rocket's - axis of symmetry as a Function of time.""" + """Aerodynamic moment in the rocket z-axis as a Function of time. + Sometimes referred to as roll moment.""" return self.__evaluate_post_process[:, [0, 12]] @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant") @@ -2225,7 +2187,7 @@ def wind_velocity_x(self): @funcify_method("Time (s)", "Wind Velocity Y (North) (m/s)", "spline", "constant") def wind_velocity_y(self): - """Wind velocity in the y direction (north) as a Function of time.""" + """Wind velocity in the Y direction (north) as a Function of time.""" return [(t, self.env.wind_velocity_y.get_value_opt(z)) for t, z in self.z] # Process fourth type of output - values calculated from previous outputs @@ -2325,17 +2287,26 @@ def path_angle(self): # Attitude Angle @funcify_method("Time (s)", "Attitude Vector X Component") def attitude_vector_x(self): - """Rocket attitude vector X component as a Function of time.""" + """Rocket attitude vector X component as a Function of time. + Same as row 1, column 3 of the rotation matrix that defines + the conversion from the body frame to the inertial frame + at each time step.""" return 2 * (self.e1 * self.e3 + self.e0 * self.e2) # a13 @funcify_method("Time (s)", "Attitude Vector Y Component") def attitude_vector_y(self): - """Rocket attitude vector Y component as a Function of time.""" + """Rocket attitude vector Y component as a Function of time. + Same as row 2, column 3 of the rotation matrix that defines + the conversion from the body frame to the inertial frame + at each time step.""" return 2 * (self.e2 * self.e3 - self.e0 * self.e1) # a23 @funcify_method("Time (s)", "Attitude Vector Z Component") def attitude_vector_z(self): - """Rocket attitude vector Z component as a Function of time.""" + """Rocket attitude vector Z component as a Function of time. + Same as row 3, column 3 of the rotation matrix that defines + the conversion from the body frame to the inertial frame + at each time step.""" return 1 - 2 * (self.e1**2 + self.e2**2) # a33 @funcify_method("Time (s)", "Attitude Angle (°)") @@ -2536,6 +2507,7 @@ def max_dynamic_pressure(self): # Total Pressure @funcify_method("Time (s)", "Total Pressure (Pa)", "spline", "zero") def total_pressure(self): + """Total pressure as a Function of time.""" return self.pressure * (1 + 0.2 * self.mach_number**2) ** (3.5) @cached_property @@ -2552,6 +2524,9 @@ def max_total_pressure(self): # Dynamics functions and variables # Aerodynamic Lift and Drag + # TODO: These are not lift and drag, they are the aerodynamic forces in + # the rocket frame, meaning they are normal and axial forces. They should + # be renamed. @funcify_method("Time (s)", "Aerodynamic Lift Force (N)", "spline", "zero") def aerodynamic_lift(self): """Aerodynamic lift force as a Function of time.""" @@ -2576,6 +2551,7 @@ def aerodynamic_spin_moment(self): # Kinetic Energy @funcify_method("Time (s)", "Rotational Kinetic Energy (J)") def rotational_energy(self): + """Rotational kinetic energy as a Function of time.""" rotational_energy = 0.5 * ( self.rocket.I_11 * self.w1**2 + self.rocket.I_22 * self.w2**2 @@ -2624,7 +2600,7 @@ def total_energy(self): # thrust Power @funcify_method("Time (s)", "thrust Power (W)", "spline", "zero") def thrust_power(self): - """thrust power as a Function of time.""" + """Thrust power as a Function of time.""" thrust = deepcopy(self.rocket.motor.thrust) thrust = thrust.set_discrete_based_on_model(self.speed) thrust_power = thrust * self.speed @@ -2639,35 +2615,105 @@ def drag_power(self): return drag_power # Angle of Attack + @cached_property + def direction_cosine_matrixes(self): + """Direction cosine matrix representing the attitude of the body frame, + relative to the inertial frame, at each time step.""" + # Stack the y_arrays from e0, e1, e2, and e3 along a new axis + stacked_arrays = np.stack( + [self.e0.y_array, self.e1.y_array, self.e2.y_array, self.e3.y_array], + axis=-1, + ) + + # Apply the transformation to the stacked arrays along the last axis + Kt = np.array([Matrix.transformation(row).transpose for row in stacked_arrays]) + + return Kt + + @cached_property + def stream_velocity_body_frame(self): + """Stream velocity array at the center of dry mass in the body frame at + each time step.""" + Kt = self.direction_cosine_matrixes + stream_velocity = np.array( + [ + self.stream_velocity_x.y_array, + self.stream_velocity_y.y_array, + self.stream_velocity_z.y_array, + ] + ).transpose() + stream_velocity_body = np.squeeze( + np.matmul(Kt, stream_velocity[:, :, np.newaxis]) + ) + return stream_velocity_body + @funcify_method("Time (s)", "Angle of Attack (°)", "spline", "constant") def angle_of_attack(self): """Angle of attack of the rocket with respect to the freestream - velocity vector.""" - dot_product = [ - -self.attitude_vector_x.get_value_opt(i) - * self.stream_velocity_x.get_value_opt(i) - - self.attitude_vector_y.get_value_opt(i) - * self.stream_velocity_y.get_value_opt(i) - - self.attitude_vector_z.get_value_opt(i) - * self.stream_velocity_z.get_value_opt(i) - for i in self.time - ] + velocity vector. Sometimes called total angle of attack. Defined as the + angle between the freestream velocity vector and the rocket's z-axis. + All in the Body Axes Coordinate System.""" + # Define stream velocity z component in body frame + stream_vz_body = ( + -self.attitude_vector_x.y_array * self.stream_velocity_x.y_array + - self.attitude_vector_y.y_array * self.stream_velocity_y.y_array + - self.attitude_vector_z.y_array * self.stream_velocity_z.y_array + ) # Define freestream speed list - free_stream_speed = [self.free_stream_speed.get_value_opt(i) for i in self.time] - free_stream_speed = np.nan_to_num(free_stream_speed) + free_stream_speed = self.free_stream_speed.y_array - # Normalize dot product - dot_product_normalized = [ - i / j if j > 1e-6 else 0 for i, j in zip(dot_product, free_stream_speed) - ] - dot_product_normalized = np.nan_to_num(dot_product_normalized) - dot_product_normalized = np.clip(dot_product_normalized, -1, 1) + stream_vz_body_normalized = np.divide( + stream_vz_body, + free_stream_speed, + out=np.zeros_like(stream_vz_body), + where=free_stream_speed > 1e-6, + ) + stream_vz_body_normalized = np.clip(stream_vz_body_normalized, -1, 1) # Calculate angle of attack and convert to degrees - angle_of_attack = np.rad2deg(np.arccos(dot_product_normalized)) + angle_of_attack = np.rad2deg(np.arccos(stream_vz_body_normalized)) + angle_of_attack = np.nan_to_num(angle_of_attack) return np.column_stack([self.time, angle_of_attack]) + @funcify_method("Time (s)", "Partial Angle of Attack (°)", "spline", "constant") + def partial_angle_of_attack(self): + """Partial angle of attack of the rocket with respect to the stream + velocity vector. By partial angle of attack, it is meant the angle + between the stream velocity vector in the y-z plane and the rocket's + z-axis. All in the Body Axes Coordinate System.""" + # Stream velocity in standard aerodynamic frame + stream_velocity = -self.stream_velocity_body_frame + alpha = np.arctan2( + stream_velocity[:, 1], + stream_velocity[:, 2], + ) # y-z plane + return np.column_stack([self.time, np.rad2deg(alpha)]) + + @funcify_method("Time (s)", "Beta (°)", "spline", "constant") + def angle_of_sideslip(self): + """Angle of sideslip of the rocket with respect to the stream + velocity vector. Defined as the angle between the stream velocity + vector in the x-z plane and the rocket's z-axis. All in the Body + Axes Coordinate System.""" + # Stream velocity in standard aerodynamic frame + stream_velocity = -self.stream_velocity_body_frame + beta = np.arctan2( + -stream_velocity[:, 0], + stream_velocity[:, 2], + ) # x-z plane + return np.column_stack([self.time, np.rad2deg(beta)]) + + @funcify_method("Time (s)", "Angle of Attack (°)", "spline", "constant") + def angle_of_attack2(self): + alpha = np.arctan( + np.sqrt( + np.tan(np.deg2rad(self.partial_angle_of_attack.y_array)) ** 2 + + np.tan(np.deg2rad(self.angle_of_sideslip.y_array)) ** 2 + ) + ) + return np.column_stack([self.time, np.rad2deg(alpha)]) + # Frequency response and stability variables @funcify_method("Frequency (Hz)", "ω1 Fourier Amplitude", "spline", "zero") def omega1_frequency_response(self): @@ -2886,9 +2932,10 @@ def __calculate_rail_button_forces(self): # TODO: complex method. # Distance from Rail Button 1 (upper) to CM rail_buttons_tuple = self.rocket.rail_buttons[0] upper_button_position = ( - rail_buttons_tuple.component.buttons_distance + rail_buttons_tuple.position + rail_buttons_tuple.component.buttons_distance + + rail_buttons_tuple.position.z ) - lower_button_position = rail_buttons_tuple.position + lower_button_position = rail_buttons_tuple.position.z angular_position_rad = rail_buttons_tuple.component.angular_position_rad D1 = ( upper_button_position - self.rocket.center_of_dry_mass_position diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 714637d1f..d9333dfc7 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -3,6 +3,7 @@ import warnings from random import choice +from rocketpy.mathutils.vector_matrix import Vector from rocketpy.motors.motor import EmptyMotor, GenericMotor, Motor from rocketpy.motors.solid_motor import SolidMotor from rocketpy.rocket.aero_surface import ( @@ -449,6 +450,9 @@ def get_surface_position(self_object, _): def _randomize_position(self, position): """Randomize a position provided as a tuple or list.""" if isinstance(position, tuple): + if isinstance(position[0], Vector): + # TODO implement randomization for X and Y positions + return position[-1](position[0].z, position[1]) return position[-1](position[0], position[1]) elif isinstance(position, list): return choice(position) if position else position @@ -496,7 +500,9 @@ def _create_surface(self, component_stochastic_surface): self.last_rnd_dict["aerodynamic_surfaces"].append( stochastic_surface.last_rnd_dict ) - self.last_rnd_dict["aerodynamic_surfaces"][-1]["position"] = position_rnd + self.last_rnd_dict["aerodynamic_surfaces"][-1]["position"] = Vector( + [0, 0, position_rnd] + ) return surface, position_rnd def _create_rail_buttons(self, component_stochastic_rail_buttons): diff --git a/tests/conftest.py b/tests/conftest.py index 6c4171b66..a23f5fd89 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -18,6 +18,8 @@ "tests.fixtures.monte_carlo.stochastic_fixtures", "tests.fixtures.monte_carlo.stochastic_motors_fixtures", "tests.fixtures.sensors.sensors_fixtures", + "tests.fixtures.generic_surfaces.generic_surfaces_fixtures", + "tests.fixtures.generic_surfaces.linear_generic_surfaces_fixtures", ] diff --git a/tests/fixtures/generic_surfaces/__init__.py b/tests/fixtures/generic_surfaces/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/fixtures/generic_surfaces/generic_surfaces_fixtures.py b/tests/fixtures/generic_surfaces/generic_surfaces_fixtures.py new file mode 100644 index 000000000..e3aa869ea --- /dev/null +++ b/tests/fixtures/generic_surfaces/generic_surfaces_fixtures.py @@ -0,0 +1,42 @@ +import pandas as pd +import pytest + + +@pytest.fixture(scope="session") +def filename_valid_coeff(tmpdir_factory): + """Creates temporary files used to test if generic surfaces + initializes correctly from CSV files""" + filename = tmpdir_factory.mktemp("aero_surface_data").join("valid_coefficients.csv") + pd.DataFrame( + { + "alpha": [0, 1, 2, 3, 0.1], + "mach": [3, 2, 1, 0, 0.2], + "cL": [4, 2, 2, 4, 5], + } + ).to_csv(filename, index=False) + + return filename + + +@pytest.fixture( + params=( + { + "alpha": [0, 1, 2, 3, 0.1], + "cL": [4, 2, 2, 4, 5], + "mach": [3, 2, 1, 0, 0.2], + }, + { + "a": [0, 1, 2, 3, 0.1], + "b": [4, 2, 2, 4, 5], + }, + ) +) +def filename_invalid_coeff(tmpdir_factory, request): + """Creates temporary CSV files used to test if generic surfaces + raises errors when initialized incorrectly from CSV files""" + filename = tmpdir_factory.mktemp("aero_surface_data").join( + "tmp_invalid_coefficients.csv" + ) + pd.DataFrame(request.param).to_csv(filename, index=False) + + return filename diff --git a/tests/fixtures/generic_surfaces/linear_generic_surfaces_fixtures.py b/tests/fixtures/generic_surfaces/linear_generic_surfaces_fixtures.py new file mode 100644 index 000000000..35ab9c1a2 --- /dev/null +++ b/tests/fixtures/generic_surfaces/linear_generic_surfaces_fixtures.py @@ -0,0 +1,44 @@ +import pandas as pd +import pytest + + +@pytest.fixture(scope="session") +def filename_valid_coeff_linear_generic_surface(tmpdir_factory): + """Creates temporary files used to test if a linear generic surface + initializes correctly from CSV files""" + filename = tmpdir_factory.mktemp("aero_surface_data").join( + "valid_coefficients_lgs.csv" + ) + pd.DataFrame( + { + "alpha": [0, 1, 2, 3, 0.1], + "mach": [3, 2, 1, 0, 0.2], + "cL_0": [4, 2, 2, 4, 5], + } + ).to_csv(filename, index=False) + + return filename + + +@pytest.fixture( + params=( + { + "alpha": [0, 1, 2, 3, 0.1], + "cL_0": [4, 2, 2, 4, 5], + "mach": [3, 2, 1, 0, 0.2], + }, + { + "a": [0, 1, 2, 3, 0.1], + "b": [4, 2, 2, 4, 5], + }, + ) +) +def filename_invalid_coeff_linear_generic_surface(tmpdir_factory, request): + """Creates temporary CSV files used to test if a linear generic surface + raises errors when initialized incorrectly from CSV files""" + filename = tmpdir_factory.mktemp("aero_surface_data").join( + "tmp_invalid_coefficients_lgs.csv" + ) + pd.DataFrame(request.param).to_csv(filename, index=False) + + return filename diff --git a/tests/fixtures/surfaces/surface_fixtures.py b/tests/fixtures/surfaces/surface_fixtures.py index 10595205f..396206bd7 100644 --- a/tests/fixtures/surfaces/surface_fixtures.py +++ b/tests/fixtures/surfaces/surface_fixtures.py @@ -1,6 +1,7 @@ import pytest from rocketpy import NoseCone, RailButtons, Tail, TrapezoidalFins +from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins @pytest.fixture @@ -62,6 +63,23 @@ def calisto_trapezoidal_fins(): ) +@pytest.fixture +def calisto_free_form_fins(): + """The free form fins of the Calisto rocket. + + Returns + ------- + rocketpy.FreeFormFins + The free form fins of the Calisto rocket. + """ + return FreeFormFins( + n=4, + shape_points=[(0, 0), (0.08, 0.1), (0.12, 0.1), (0.12, 0)], + rocket_radius=0.0635, + name="calisto_free_form_fins", + ) + + @pytest.fixture def calisto_rail_buttons(): """The rail buttons of the Calisto rocket. diff --git a/tests/integration/test_flight.py b/tests/integration/test_flight.py index 422941c4d..fc1dd1956 100644 --- a/tests/integration/test_flight.py +++ b/tests/integration/test_flight.py @@ -335,8 +335,8 @@ def test_rolling_flight( # pylint: disable=unused-argument position=-1.04956, cant_angle=0.5, ) - calisto.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) - calisto.aerodynamic_surfaces.add(calisto_tail, -1.313) + calisto.add_surfaces(calisto_nose_cone, 1.160) + calisto.add_surfaces(calisto_tail, -1.313) calisto.parachutes.append(calisto_main_chute) calisto.parachutes.append(calisto_drogue_chute) @@ -365,9 +365,9 @@ def test_eccentricity_on_flight( # pylint: disable=unused-argument test_rocket.set_rail_buttons(0.082, -0.618) test_rocket.add_motor(cesaroni_m1670, position=-1.373) - calisto.aerodynamic_surfaces.add(calisto_trapezoidal_fins, position=-1.04956) - calisto.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) - calisto.aerodynamic_surfaces.add(calisto_tail, -1.313) + calisto.add_surfaces(calisto_trapezoidal_fins, -1.04956) + calisto.add_surfaces(calisto_nose_cone, 1.160) + calisto.add_surfaces(calisto_tail, -1.313) calisto.add_cm_eccentricity(x=-0.01, y=-0.01) test_flight = Flight( diff --git a/tests/integration/test_rocket.py b/tests/integration/test_rocket.py index 4d5daf7a6..e51fe925e 100644 --- a/tests/integration/test_rocket.py +++ b/tests/integration/test_rocket.py @@ -1,6 +1,7 @@ from unittest.mock import patch import numpy as np +import pytest @patch("matplotlib.pyplot.show") @@ -14,8 +15,8 @@ def test_airfoil( ): test_rocket = calisto test_rocket.set_rail_buttons(0.082, -0.618) - calisto.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) - calisto.aerodynamic_surfaces.add(calisto_tail, -1.313) + calisto.add_surfaces(calisto_nose_cone, 1.160) + calisto.add_surfaces(calisto_tail, -1.313) test_rocket.add_trapezoidal_fins( 2, @@ -126,9 +127,47 @@ def test_rocket(mock_show, calisto_robust): # pylint: disable=unused-argument @patch("matplotlib.pyplot.show") def test_aero_surfaces_infos( # pylint: disable=unused-argument - mock_show, calisto_nose_cone, calisto_tail, calisto_trapezoidal_fins + mock_show, + calisto_nose_cone, + calisto_tail, + calisto_trapezoidal_fins, + calisto_free_form_fins, ): assert calisto_nose_cone.all_info() is None assert calisto_trapezoidal_fins.all_info() is None assert calisto_tail.all_info() is None - assert calisto_trapezoidal_fins.draw() is None + assert calisto_trapezoidal_fins.all_info() is None + assert calisto_free_form_fins.all_info() is None + + +@pytest.mark.parametrize( + "attribute, tolerance", + [ + ("cpz", 1e-3), + ("clalpha", 1e-3), + ("Af", 1e-3), + ("gamma_c", 1e-3), + ("Yma", 1e-3), + ("tau", 1e-3), + ("roll_geometrical_constant", 1e-3), + ("lift_interference_factor", 1e-3), + ("roll_forcing_interference_factor", 1e-3), + ("roll_damping_interference_factor", 1e-2), + ], +) +def test_calisto_free_form_fins_equivalence( + calisto_free_form_fins, calisto_trapezoidal_fins, attribute, tolerance +): + """Test the equivalence of the free form fins with the same geometric + characteristics as the trapezoidal fins, comparing cp, cnalpha, and + geometrical parameters.""" + + # Handle the 'clalpha' method comparison differently as it's a callable + if attribute == "clalpha": + free_form_value = calisto_free_form_fins.clalpha(0) + trapezoidal_value = calisto_trapezoidal_fins.clalpha(0) + else: + free_form_value = getattr(calisto_free_form_fins, attribute) + trapezoidal_value = getattr(calisto_trapezoidal_fins, attribute) + + assert abs(free_form_value - trapezoidal_value) < tolerance diff --git a/tests/unit/test_environment.py b/tests/unit/test_environment.py index 7769f7b85..714470a40 100644 --- a/tests/unit/test_environment.py +++ b/tests/unit/test_environment.py @@ -71,13 +71,18 @@ def test_location_set_topographic_profile_computes_elevation( def test_geodesic_coordinate_geodesic_to_utm_converts_coordinate(): """Tests the conversion from geodesic to UTM coordinates.""" - x, y, utm_zone, utm_letter, north_south_hemis, east_west_hemis = ( - Environment.geodesic_to_utm( - lat=32.990254, - lon=-106.974998, - semi_major_axis=6378137.0, # WGS84 - flattening=1 / 298.257223563, # WGS84 - ) + ( + x, + y, + utm_zone, + utm_letter, + north_south_hemis, + east_west_hemis, + ) = Environment.geodesic_to_utm( + lat=32.990254, + lon=-106.974998, + semi_major_axis=6378137.0, # WGS84 + flattening=1 / 298.257223563, # WGS84 ) assert np.isclose(x, 315468.64, atol=1e-5) assert np.isclose(y, 3651938.65, atol=1e-5) diff --git a/tests/unit/test_generic_surfaces.py b/tests/unit/test_generic_surfaces.py new file mode 100644 index 000000000..a04ed429d --- /dev/null +++ b/tests/unit/test_generic_surfaces.py @@ -0,0 +1,91 @@ +import pytest + +from rocketpy import Function, GenericSurface +from rocketpy.mathutils import Vector + +REFERENCE_AREA = 1 +REFERENCE_LENGTH = 1 + + +@pytest.mark.parametrize( + "coefficients", + [ + "cL", + {"invalid_name": 0}, + {"cL": "inexistent_file.csv"}, + {"cL": Function(lambda x1, x2, x3, x4, x5, x6: 0)}, + {"cL": lambda x1: 0}, + {"cL": {}}, + ], +) +def test_invalid_initialization(coefficients): + """Checks if generic surface raises errors in initialization + when coefficient argument is invalid""" + + with pytest.raises((ValueError, TypeError)): + GenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients=coefficients, + ) + + +def test_invalid_initialization_from_csv(filename_invalid_coeff): + """Checks if generic surfaces raises errors when initialized incorrectly + from a csv file""" + with pytest.raises(ValueError): + GenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients={"cL": str(filename_invalid_coeff)}, + ) + + +@pytest.mark.parametrize( + "coefficients", + [ + {}, + {"cL": 0}, + { + "cL": 0, + "cQ": Function(lambda x1, x2, x3, x4, x5, x6, x7: 0), + "cD": lambda x1, x2, x3, x4, x5, x6, x7: 0, + }, + ], +) +def test_valid_initialization(coefficients): + """Checks if generic surface initializes correctly when coefficient + argument is valid""" + + GenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients=coefficients, + ) + + +def test_valid_initialization_from_csv(filename_valid_coeff): + """Checks if generic surfaces initializes correctly when + coefficients is set from a csv file""" + GenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients={"cL": str(filename_valid_coeff)}, + ) + + +def test_compute_forces_and_moments(): + """Checks if there are not logical errors in + compute forces and moments""" + + gs_object = GenericSurface(REFERENCE_AREA, REFERENCE_LENGTH, {}) + forces_and_moments = gs_object.compute_forces_and_moments( + stream_velocity=Vector((0, 0, 0)), + stream_speed=0, + stream_mach=0, + rho=0, + cp=Vector((0, 0, 0)), + omega=(0, 0, 0), + reynolds=0, + ) + assert forces_and_moments == (0, 0, 0, 0, 0, 0) diff --git a/tests/unit/test_linear_generic_surfaces.py b/tests/unit/test_linear_generic_surfaces.py new file mode 100644 index 000000000..858226444 --- /dev/null +++ b/tests/unit/test_linear_generic_surfaces.py @@ -0,0 +1,91 @@ +import pytest + +from rocketpy import Function, LinearGenericSurface +from rocketpy.mathutils import Vector + +REFERENCE_AREA = 1 +REFERENCE_LENGTH = 1 + + +@pytest.mark.parametrize( + "coefficients", + [ + "cL_0", + {"invalid_name": 0}, + {"cL_0": "inexistent_file.csv"}, + {"cL_0": Function(lambda x1, x2, x3, x4, x5, x6: 0)}, + {"cL_0": lambda x1: 0}, + {"cL_0": {}}, + ], +) +def test_invalid_initialization(coefficients): + """Checks if linear generic surface raises errors in initialization + when coefficient argument is invalid""" + + with pytest.raises((ValueError, TypeError)): + LinearGenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients=coefficients, + ) + + +def test_invalid_initialization_from_csv(filename_invalid_coeff_linear_generic_surface): + """Checks if linear generic surfaces raises errors when initialized incorrectly + from a csv file""" + with pytest.raises(ValueError): + LinearGenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients={"cL_0": str(filename_invalid_coeff_linear_generic_surface)}, + ) + + +@pytest.mark.parametrize( + "coefficients", + [ + {}, + {"cL_0": 0}, + { + "cL_0": 0, + "cQ_0": Function(lambda x1, x2, x3, x4, x5, x6, x7: 0), + "cD_0": lambda x1, x2, x3, x4, x5, x6, x7: 0, + }, + ], +) +def test_valid_initialization(coefficients): + """Checks if linear generic surface initializes correctly when coefficient + argument is valid""" + + LinearGenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients=coefficients, + ) + + +def test_valid_initialization_from_csv(filename_valid_coeff_linear_generic_surface): + """Checks if linear generic surfaces initializes correctly when + coefficients is set from a csv file""" + LinearGenericSurface( + reference_area=REFERENCE_AREA, + reference_length=REFERENCE_LENGTH, + coefficients={"cL_0": str(filename_valid_coeff_linear_generic_surface)}, + ) + + +def test_compute_forces_and_moments(): + """Checks if there are not logical errors in + compute forces and moments""" + + lgs_object = LinearGenericSurface(REFERENCE_AREA, REFERENCE_LENGTH, {}) + forces_and_moments = lgs_object.compute_forces_and_moments( + stream_velocity=Vector((0, 0, 0)), + stream_speed=1, + stream_mach=0, + rho=0, + cp=Vector((0, 0, 0)), + omega=(0, 0, 0), + reynolds=0, + ) + assert forces_and_moments == (0, 0, 0, 0, 0, 0) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index a0d9b79a3..f048f4b3c 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -4,6 +4,7 @@ import pytest from rocketpy import Function, NoseCone, Rocket, SolidMotor +from rocketpy.mathutils.vector_matrix import Vector from rocketpy.motors.motor import EmptyMotor, Motor @@ -124,7 +125,7 @@ def test_add_trapezoidal_fins_sweep_angle( calisto_nose_cone, ): # Reference values from OpenRocket - calisto.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) + calisto.add_surfaces(calisto_nose_cone, Vector([0, 0, 1.160])) fin_set = calisto.add_trapezoidal_fins( n=3, span=0.090, @@ -164,7 +165,7 @@ def test_add_trapezoidal_fins_sweep_length( calisto_nose_cone, ): # Reference values from OpenRocket - calisto.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) + calisto.add_surfaces(calisto_nose_cone, Vector([0, 0, 1.160])) fin_set = calisto.add_trapezoidal_fins( n=3, span=0.090, @@ -379,7 +380,7 @@ def test_set_rail_button(calisto): == pytest.approx(0.7, 1e-12) ) # assert buttons position on rocket - assert calisto.rail_buttons[0].position == -0.5 + assert calisto.rail_buttons[0].position.z == -0.5 # assert angular position assert ( rail_buttons.angular_position @@ -389,16 +390,7 @@ def test_set_rail_button(calisto): # assert upper button position assert calisto.rail_buttons[0].component.buttons_distance + calisto.rail_buttons[ 0 - ].position == pytest.approx(0.2, 1e-12) - - -def test_add_rail_button(calisto, calisto_rail_buttons): - calisto.add_surfaces(calisto_rail_buttons, -0.5) - assert calisto.rail_buttons[0].position == -0.5 - upper_position = ( - calisto_rail_buttons.buttons_distance + calisto.rail_buttons[0].position - ) - assert upper_position == pytest.approx(0.2, 1e-12) + ].position.z == pytest.approx(0.2, 1e-12) def test_evaluate_total_mass(calisto_motorless): @@ -641,8 +633,8 @@ def test_coordinate_system_orientation( rocket_tail_to_nose.add_motor(motor_nozzle_to_combustion_chamber, position=-1.373) - rocket_tail_to_nose.aerodynamic_surfaces.add(calisto_nose_cone, 1.160) - rocket_tail_to_nose.aerodynamic_surfaces.add(calisto_trapezoidal_fins, -1.168) + rocket_tail_to_nose.add_surfaces(calisto_nose_cone, 1.160) + rocket_tail_to_nose.add_surfaces(calisto_trapezoidal_fins, -1.168) static_margin_tail_to_nose = rocket_tail_to_nose.static_margin @@ -658,8 +650,8 @@ def test_coordinate_system_orientation( rocket_nose_to_tail.add_motor(motor_combustion_chamber_to_nozzle, position=1.373) - rocket_nose_to_tail.aerodynamic_surfaces.add(calisto_nose_cone, -1.160) - rocket_nose_to_tail.aerodynamic_surfaces.add(calisto_trapezoidal_fins, 1.168) + rocket_nose_to_tail.add_surfaces(calisto_nose_cone, -1.160) + rocket_nose_to_tail.add_surfaces(calisto_trapezoidal_fins, 1.168) static_margin_nose_to_tail = rocket_nose_to_tail.static_margin