ENH: Introduce Net Thrust with pressure corrections

rocketpy/motors/empty_motor.py

@@ -13,6 +13,7 @@ def __init__(self):
             dry_inertia=(0, 0, 0),
             nozzle_radius=0,
             center_of_dry_mass_position=0,
+            reference_pressure=0,
             dry_mass=0,
             nozzle_position=0,
             burn_time=1,

rocketpy/motors/hybrid_motor.py

@@ -29,6 +29,8 @@ class HybridMotor(Motor):
         "combustion_chamber_to_nozzle".
     HybridMotor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    HybridMotor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     HybridMotor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See
@@ -147,7 +149,11 @@ class HybridMotor(Motor):
         HybridMotor.propellant_I_22 and HybridMotor.propellant_I_33 for
         more information.
     HybridMotor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from thrust source, in Newtons, as a
+        function of time.
+    HybridMotor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     HybridMotor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     HybridMotor.max_thrust : float
@@ -181,6 +187,9 @@ class HybridMotor(Motor):
         Method of interpolation used in case thrust curve is given
         by data set in .csv or .eng, or as an array. Options are 'spline'
         'akima' and 'linear'. Default is "linear".
+    HybridMotor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     def __init__(  # pylint: disable=too-many-arguments
@@ -197,6 +206,7 @@ def __init__(  # pylint: disable=too-many-arguments
         grain_separation,
         grains_center_of_mass_position,
         center_of_dry_mass_position,
+        reference_pressure=None,
         nozzle_position=0,
         burn_time=None,
         throat_radius=0.01,
@@ -258,6 +268,8 @@ class Function. Thrust units are Newtons.
             The position, in meters, of the motor's center of mass with respect
             to the motor's coordinate system when it is devoid of propellant.
             See :doc:`Positions and Coordinate Systems </user/positions>`.
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         nozzle_position : int, float, optional
             Motor's nozzle outlet position in meters, in the motor's coordinate
             system. See :doc:`Positions and Coordinate Systems </user/positions>`
@@ -308,6 +320,7 @@ class Function. Thrust units are Newtons.
             dry_inertia=dry_inertia,
             nozzle_radius=nozzle_radius,
             center_of_dry_mass_position=center_of_dry_mass_position,
+            reference_pressure=reference_pressure,
             dry_mass=dry_mass,
             nozzle_position=nozzle_position,
             burn_time=burn_time,
@@ -321,6 +334,7 @@ class Function. Thrust units are Newtons.
             dry_inertia,
             nozzle_radius,
             center_of_dry_mass_position,
+            reference_pressure,
             nozzle_position,
             burn_time,
             reshape_thrust_curve,
@@ -340,6 +354,7 @@ class Function. Thrust units are Newtons.
             grain_separation,
             grains_center_of_mass_position,
             center_of_dry_mass_position,
+            reference_pressure,
             nozzle_position,
             burn_time,
             throat_radius,
@@ -657,6 +672,7 @@ def from_dict(cls, data):
             nozzle_radius=data["nozzle_radius"],
             dry_mass=data["dry_mass"],
             center_of_dry_mass_position=data["center_of_dry_mass_position"],
+            reference_pressure=data["reference_pressure"],
             dry_inertia=(
                 data["dry_I_11"],
                 data["dry_I_22"],

rocketpy/motors/liquid_motor.py

@@ -29,6 +29,8 @@ class LiquidMotor(Motor):
         "combustion_chamber_to_nozzle".
     LiquidMotor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    LiquidMotor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     LiquidMotor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See
@@ -122,7 +124,11 @@ class LiquidMotor(Motor):
         LiquidMotor.propellant_I_22 and LiquidMotor.propellant_I_33 for
         more information.
     LiquidMotor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from thrust source, in Newtons, as a
+        function of time.
+    LiquidMotor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     LiquidMotor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     LiquidMotor.max_thrust : float
@@ -143,6 +149,9 @@ class LiquidMotor(Motor):
         burn_out_time and the burn_start_time.
     LiquidMotor.exhaust_velocity : Function
         Propulsion gases exhaust velocity in m/s.
+    LiquidMotor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     def __init__(
@@ -152,6 +161,7 @@ def __init__(
         dry_inertia,
         nozzle_radius,
         center_of_dry_mass_position,
+        reference_pressure=None,
         nozzle_position=0,
         burn_time=None,
         reshape_thrust_curve=False,
@@ -194,6 +204,8 @@ class Function. Thrust units are Newtons.
             The position, in meters, of the motor's center of mass with respect
             to the motor's coordinate system when it is devoid of propellant.
             See :doc:`Positions and Coordinate Systems </user/positions>`
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         nozzle_position : float
             Motor's nozzle outlet position in meters, specified in the motor's
             coordinate system. See
@@ -236,6 +248,7 @@ class Function. Thrust units are Newtons.
             dry_inertia=dry_inertia,
             nozzle_radius=nozzle_radius,
             center_of_dry_mass_position=center_of_dry_mass_position,
+            reference_pressure=reference_pressure,
             dry_mass=dry_mass,
             nozzle_position=nozzle_position,
             burn_time=burn_time,
@@ -500,6 +513,7 @@ def from_dict(cls, data):
             nozzle_radius=data["nozzle_radius"],
             dry_mass=data["dry_mass"],
             center_of_dry_mass_position=data["center_of_dry_mass_position"],
+            reference_pressure=data["reference_pressure"],
             dry_inertia=(
                 data["dry_I_11"],
                 data["dry_I_22"],

rocketpy/motors/motor.py

@@ -28,6 +28,8 @@ class Motor(ABC):
         "combustion_chamber_to_nozzle".
     Motor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    Motor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     Motor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See :ref:`positions` for more information.
@@ -121,7 +123,11 @@ class Motor(ABC):
         e_3 axes in kg*m^2, as a function of time. See Motor.propellant_I_22
         and Motor.propellant_I_33 for more information.
     Motor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from the thrust source, in Newtons, as a
+        function of time.
+    Motor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     Motor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     Motor.max_thrust : float
@@ -146,6 +152,9 @@ class Motor(ABC):
         Method of interpolation used in case thrust curve is given
         by data set in .csv or .eng, or as an array. Options are 'spline'
         'akima' and 'linear'. Default is "linear".
+    Motor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     # pylint: disable=too-many-statements
@@ -155,6 +164,7 @@ def __init__(
         dry_inertia,
         nozzle_radius,
         center_of_dry_mass_position,
+        reference_pressure=None,
         dry_mass=None,
         nozzle_position=0,
         burn_time=None,
@@ -199,6 +209,8 @@ class Function. Thrust units are Newtons.
             (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0.
         nozzle_radius : int, float, optional
             Motor's nozzle outlet radius in meters.
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         burn_time: float, tuple of float, optional
             Motor's burn time.
             If a float is given, the burn time is assumed to be between 0 and
@@ -257,7 +269,9 @@ class Function. Thrust units are Newtons.
         self.interpolate = interpolation_method
         self.nozzle_position = nozzle_position
         self.nozzle_radius = nozzle_radius
+        self.nozzle_area = np.pi * nozzle_radius**2
         self.center_of_dry_mass_position = center_of_dry_mass_position
+        self.reference_pressure = reference_pressure
 
         # Inertia tensor setup
         inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia
@@ -302,15 +316,27 @@ class Function. Thrust units are Newtons.
         # Post process thrust
         self.thrust = Motor.clip_thrust(self.thrust, self.burn_time)
 
+        # Evaluate vacuum thrust
+        if reference_pressure is None:
+            vacuum_thrust_source = self.thrust.source
+        else:
+            vacuum_thrust_source = self.get_vacuum_thrust(
+                self.thrust.source, reference_pressure
+            )
+        self.vacuum_thrust_source = vacuum_thrust_source
+        self.vacuum_thrust = Function(
+            vacuum_thrust_source,
+            "Time (s)",
+            "Vacuum Thrust (N)",
+            self.interpolate,
+            "zero",
+        )
+
         # Auxiliary quantities
         self.burn_start_time = self.burn_time[0]
         self.burn_out_time = self.burn_time[1]
         self.burn_duration = self.burn_time[1] - self.burn_time[0]
 
-        # Define motor attributes
-        self.nozzle_radius = nozzle_radius
-        self.nozzle_position = nozzle_position
-
         # Compute thrust metrics
         self.max_thrust = np.amax(self.thrust.y_array)
         max_thrust_index = np.argmax(self.thrust.y_array)
@@ -1039,6 +1065,34 @@ def import_eng(file_name):
         # Return all extract content
         return comments, description, data_points
 
+    def get_vacuum_thrust(self, thrust_source, reference_pressure):
+        """Calculate the vacuum thrust from the raw thrust and the reference
+        pressure at which the thrust data was recorded.
+
+        Parameters
+        ----------
+        thrust_source : list
+            A list of points representing the thrust curve.
+        reference_pressure : int, float
+            The atmospheric pressure at which the thrust data was recorded.
+
+        Returns
+        -------
+        vacuum_thrust_source : list
+            A list of points representing the vacuum thrust curve.
+        """
+        # Initialize arrays
+        vacuum_thrust_source = []
+
+        # Reduction in thrust due to atmospheric pressure
+        thrust_reduction = reference_pressure * self.nozzle_area
+        for point in thrust_source:
+            time = point[0]
+            thrust = point[1]
+            vacuum_thrust_source.append([time, thrust + thrust_reduction])
+
+        return vacuum_thrust_source
+
     def export_eng(self, file_name, motor_name):
         """Exports thrust curve data points and motor description to
         .eng file format. A description of the format can be found
@@ -1105,6 +1159,8 @@ def to_dict(self, include_outputs=False):
             "dry_I_13": self.dry_I_13,
             "dry_I_23": self.dry_I_23,
             "nozzle_radius": self.nozzle_radius,
+            "nozzle_area": self.nozzle_area,
+            "reference_pressure": self.reference_pressure,
             "center_of_dry_mass_position": self.center_of_dry_mass_position,
             "dry_mass": self.dry_mass,
             "nozzle_position": self.nozzle_position,
@@ -1116,6 +1172,7 @@ def to_dict(self, include_outputs=False):
         if include_outputs:
             data.update(
                 {
+                    "vacuum_thrust": self.vacuum_thrust,
                     "total_mass": self.total_mass,
                     "propellant_mass": self.propellant_mass,
                     "mass_flow_rate": self.mass_flow_rate,
@@ -1176,6 +1233,7 @@ class GenericMotor(Motor):
     therefore for more accurate results, use the ``SolidMotor``, ``HybridMotor``
     or ``LiquidMotor`` classes."""
 
+    # pylint: disable=too-many-arguments
     def __init__(
         self,
         thrust_source,
@@ -1185,6 +1243,7 @@ def __init__(
         chamber_position,
         propellant_initial_mass,
         nozzle_radius,
+        reference_pressure=None,
         dry_mass=0,
         center_of_dry_mass_position=None,
         dry_inertia=(0, 0, 0),
@@ -1231,6 +1290,8 @@ def __init__(
             to the motor's coordinate system when it is devoid of propellant.
             If not specified, automatically sourced as the chamber position.
             See :doc:`Positions and Coordinate Systems </user/positions>`
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         dry_inertia : tuple, list
             Tuple or list containing the motor's dry mass inertia tensor
             components, in kg*m^2. This inertia is defined with respect to the
@@ -1287,6 +1348,7 @@ def __init__(
             dry_inertia=dry_inertia,
             nozzle_radius=nozzle_radius,
             center_of_dry_mass_position=center_of_dry_mass_position,
+            reference_pressure=reference_pressure,
             dry_mass=dry_mass,
             nozzle_position=nozzle_position,
             burn_time=burn_time,
@@ -1449,6 +1511,7 @@ def load_from_eng_file(
         dry_mass=None,
         burn_time=None,
         center_of_dry_mass_position=None,
+        reference_pressure=None,
         dry_inertia=(0, 0, 0),
         nozzle_position=0,
         reshape_thrust_curve=False,
@@ -1490,6 +1553,8 @@ def load_from_eng_file(
             The position, in meters, of the motor's center of mass with respect
             to the motor's coordinate system when it is devoid of propellant.
             If not specified, automatically sourced as the chamber position.
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         dry_inertia : tuple, list
             Tuple or list containing the motor's dry mass inertia tensor
         nozzle_position : int, float, optional
@@ -1557,6 +1622,7 @@ def load_from_eng_file(
             chamber_position=chamber_position,
             propellant_initial_mass=propellant_initial_mass,
             nozzle_radius=nozzle_radius,
+            reference_pressure=reference_pressure,
             dry_mass=dry_mass,
             center_of_dry_mass_position=center_of_dry_mass_position,
             dry_inertia=dry_inertia,
@@ -1594,6 +1660,7 @@ def from_dict(cls, data):
             chamber_position=data["chamber_position"],
             propellant_initial_mass=data["propellant_initial_mass"],
             nozzle_radius=data["nozzle_radius"],
+            reference_pressure=data["reference_pressure"],
             dry_mass=data["dry_mass"],
             center_of_dry_mass_position=data["center_of_dry_mass_position"],
             dry_inertia=(