ENH: Implement Parachute Opening Shock Force Estimation

rocketpy/rocket/parachute.py

@@ -119,6 +119,7 @@ def __init__(
         radius=1.5,
         height=None,
         porosity=0.0432,
+        opening_shock_coefficient=1.6,  # TODO: analyse methods to calculate Cx and X1
     ):
         """Initializes Parachute class.
 
@@ -184,6 +185,14 @@ def __init__(
             in [0, 1]. Affects only the added-mass scaling during descent; it does
             not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass
             of 1.0 (“neutral” behavior).
+        opening_shock_coefficient : float, optional
+            Coefficient used to calculate the parachute's opening shock force.
+            Combined factor (Cx * X1) accounting for canopy shape and mass ratio.
+            Default value set for 1.6, can be calculated by geometrical and porosity
+            relations.
+        opening_shock_force : float, optional
+            Parachute's estimated opening shock force.
+            Calculated based on Knacke, T. W. (1992). Parachute Recovery Systems Design Manual.
         """
         self.name = name
         self.cd_s = cd_s
@@ -203,6 +212,8 @@ def __init__(
         self.radius = radius
         self.height = height or radius
         self.porosity = porosity
+        self.opening_shock_coefficient = opening_shock_coefficient
+        self.opening_shock_force = None
         self.added_mass_coefficient = 1.068 * (
             1
             - 1.465 * self.porosity
@@ -310,6 +321,7 @@ def to_dict(self, **kwargs):
             "radius": self.radius,
             "height": self.height,
             "porosity": self.porosity,
+            "opening_shock_coefficient": self.opening_shock_coefficient,
         }
 
         if kwargs.get("include_outputs", False):
@@ -321,6 +333,7 @@ def to_dict(self, **kwargs):
             )
             data["noisy_pressure_signal"] = self.noisy_pressure_signal
             data["clean_pressure_signal"] = self.clean_pressure_signal
+            data["opening_shock_force"] = self.opening_shock_force
 
         return data
 
@@ -343,6 +356,32 @@ def from_dict(cls, data):
             radius=data.get("radius", 1.5),
             height=data.get("height", None),
             porosity=data.get("porosity", 0.0432),
+            opening_shock_coefficient=data.get("opening_shock_coefficient", 1.6),
         )
 
         return parachute
+
+    def calculate_opening_shock(self, density, velocity):
+        """
+        Calculates the opening shock force based on Knacke's formula.
+
+        Fo = Cx * X1 * q * S * Cd
+        Knacke, T. W. (1992). Parachute Recovery Systems Design Manual.(Page 5-50)
+
+        Parameters
+        ----------
+        density: float
+            Air density during the parachute's opening (kg/m^3).
+        velocity: float
+            Rocket velocity relative to the air during the parachute's opening (m/s).
+
+        Returns
+        -------
+        force: float
+            The estimated peak opening shock force during the parachute's opening (N).
+        """
+
+        dynamic_pressure = 0.5 * density * (velocity**2)
+        force = self.opening_shock_coefficient * dynamic_pressure * self.cd_s
+
+        return force

rocketpy/simulation/flight.py

@@ -773,6 +773,10 @@ def __simulate(self, verbose):
                         self.y_sol,
                         self.sensors,
                     ):
+
+                        # Calculates the parachute's opening shock force
+                        self.calculate_parachute_opening_shock_force(parachute)
+
                         # Remove parachute from flight parachutes
                         self.parachutes.remove(parachute)
                         # Create phase for time after detection and before inflation
@@ -800,8 +804,7 @@ def __simulate(self, verbose):
                             lambda self, parachute_porosity=parachute.porosity: setattr(
                                 self, "parachute_porosity", parachute_porosity
                             ),
-                            lambda self,
-                            added_mass_coefficient=parachute.added_mass_coefficient: setattr(
+                            lambda self, added_mass_coefficient=parachute.added_mass_coefficient: setattr(
                                 self,
                                 "parachute_added_mass_coefficient",
                                 added_mass_coefficient,
@@ -1048,30 +1051,25 @@ def __simulate(self, verbose):
                                             i += 1
                                         # Create flight phase for time after inflation
                                         callbacks = [
-                                            lambda self,
-                                            parachute_cd_s=parachute.cd_s: setattr(
+                                            lambda self, parachute_cd_s=parachute.cd_s: setattr(
                                                 self, "parachute_cd_s", parachute_cd_s
                                             ),
-                                            lambda self,
-                                            parachute_radius=parachute.radius: setattr(
+                                            lambda self, parachute_radius=parachute.radius: setattr(
                                                 self,
                                                 "parachute_radius",
                                                 parachute_radius,
                                             ),
-                                            lambda self,
-                                            parachute_height=parachute.height: setattr(
+                                            lambda self, parachute_height=parachute.height: setattr(
                                                 self,
                                                 "parachute_height",
                                                 parachute_height,
                                             ),
-                                            lambda self,
-                                            parachute_porosity=parachute.porosity: setattr(
+                                            lambda self, parachute_porosity=parachute.porosity: setattr(
                                                 self,
                                                 "parachute_porosity",
                                                 parachute_porosity,
                                             ),
-                                            lambda self,
-                                            added_mass_coefficient=parachute.added_mass_coefficient: setattr(
+                                            lambda self, added_mass_coefficient=parachute.added_mass_coefficient: setattr(
                                                 self,
                                                 "parachute_added_mass_coefficient",
                                                 added_mass_coefficient,
@@ -1098,6 +1096,12 @@ def __simulate(self, verbose):
                                         phase.time_nodes.flush_after(node_index)
                                         phase.time_nodes.add_node(self.t, [], [], [])
                                         phase.solver.status = "finished"
+
+                                        # Calculates the parachute's opening shock force
+                                        self.calculate_parachute_opening_shock_force(
+                                            parachute
+                                        )
+
                                         # Save parachute event
                                         self.parachute_events.append(
                                             [self.t, parachute]
@@ -1601,7 +1605,9 @@ def udot_rail2(self, t, u, post_processing=False):  # pragma: no cover
         # Hey! We will finish this function later, now we just can use u_dot
         return self.u_dot_generalized(t, u, post_processing=post_processing)
 
-    def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals,too-many-statements
+    def u_dot(
+        self, t, u, post_processing=False
+    ):  # pylint: disable=too-many-locals,too-many-statements
         """Calculates derivative of u state vector with respect to time
         when rocket is flying in 6 DOF motion during ascent out of rail
         and descent without parachute.
@@ -2147,7 +2153,9 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False):
 
         return u_dot
 
-    def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too-many-locals,too-many-statements
+    def u_dot_generalized(
+        self, t, u, post_processing=False
+    ):  # pylint: disable=too-many-locals,too-many-statements
         """Calculates derivative of u state vector with respect to time when the
         rocket is flying in 6 DOF motion in space and significant mass variation
         effects exist. Typical flight phases include powered ascent after launch
@@ -3982,9 +3990,7 @@ def add(self, flight_phase, index=None):  # TODO: quite complex method
             new_index = (
                 index - 1
                 if flight_phase.t < previous_phase.t
-                else index + 1
-                if flight_phase.t > next_phase.t
-                else index
+                else index + 1 if flight_phase.t > next_phase.t else index
             )
             flight_phase.t += adjust
             self.add(flight_phase, new_index)
@@ -4359,3 +4365,18 @@ def rail_button2_bending_moment(self):
     def max_rail_button2_bending_moment(self):
         """Maximum lower rail button bending moment, in N·m."""
         return self.calculate_rail_button_bending_moments[3]
+
+    def calculate_parachute_opening_shock_force(self, parachute):
+        """Calculates and stores the shock force on parachute opening
+        Uses the current self.y_sol and self.env.
+        """
+        # Calculate opening shock force
+        opening_altitude = self.y_sol[2]
+        opening_density = self.env.density(opening_altitude)
+        opening_velocity = (
+            (self.y_sol[3]) ** 2 + (self.y_sol[4]) ** 2 + (self.y_sol[5]) ** 2
+        ) ** 0.5
+
+        parachute.opening_shock_force = parachute.calculate_opening_shock(
+            opening_density, opening_velocity
+        )

tests/integration/simulation/test_flight.py

@@ -4,7 +4,7 @@
 import numpy as np
 import pytest
 
-from rocketpy import Flight
+from rocketpy import Flight, Parachute
 
 plt.rcParams.update({"figure.max_open_warning": 0})
 
@@ -69,7 +69,9 @@ def test_all_info_different_solvers(
 
 
 @patch("matplotlib.pyplot.show")
-def test_hybrid_motor_flight(mock_show, flight_calisto_hybrid_modded):  # pylint: disable=unused-argument
+def test_hybrid_motor_flight(
+    mock_show, flight_calisto_hybrid_modded
+):  # pylint: disable=unused-argument
     """Test the flight of a rocket with a hybrid motor. This test only validates
     that a flight simulation can be performed with a hybrid motor; it does not
     validate the results.
@@ -85,7 +87,9 @@ def test_hybrid_motor_flight(mock_show, flight_calisto_hybrid_modded):  # pylint
 
 
 @patch("matplotlib.pyplot.show")
-def test_liquid_motor_flight(mock_show, flight_calisto_liquid_modded):  # pylint: disable=unused-argument
+def test_liquid_motor_flight(
+    mock_show, flight_calisto_liquid_modded
+):  # pylint: disable=unused-argument
     """Test the flight of a rocket with a liquid motor. This test only validates
     that a flight simulation can be performed with a liquid motor; it does not
     validate the results.
@@ -102,7 +106,9 @@ def test_liquid_motor_flight(mock_show, flight_calisto_liquid_modded):  # pylint
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_time_overshoot(mock_show, calisto_robust, example_spaceport_env):  # pylint: disable=unused-argument
+def test_time_overshoot(
+    mock_show, calisto_robust, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Test the time_overshoot parameter of the Flight class. This basically
     calls the all_info() method for a simulation without time_overshoot and
     checks if it returns None. It is not testing if the values are correct,
@@ -131,7 +137,9 @@ def test_time_overshoot(mock_show, calisto_robust, example_spaceport_env):  # py
 
 
 @patch("matplotlib.pyplot.show")
-def test_simpler_parachute_triggers(mock_show, example_plain_env, calisto_robust):  # pylint: disable=unused-argument
+def test_simpler_parachute_triggers(
+    mock_show, example_plain_env, calisto_robust
+):  # pylint: disable=unused-argument
     """Tests different types of parachute triggers. This is important to ensure
     the code is working as intended, since the parachute triggers can have very
     different format definitions. It will add 3 parachutes using different
@@ -273,7 +281,9 @@ def test_eccentricity_on_flight(  # pylint: disable=unused-argument
 
 
 @patch("matplotlib.pyplot.show")
-def test_air_brakes_flight(mock_show, flight_calisto_air_brakes):  # pylint: disable=unused-argument
+def test_air_brakes_flight(
+    mock_show, flight_calisto_air_brakes
+):  # pylint: disable=unused-argument
     """Test the flight of a rocket with air brakes. This test only validates
     that a flight simulation can be performed with air brakes; it does not
     validate the results.
@@ -294,7 +304,9 @@ def test_air_brakes_flight(mock_show, flight_calisto_air_brakes):  # pylint: dis
 
 
 @patch("matplotlib.pyplot.show")
-def test_initial_solution(mock_show, example_plain_env, calisto_robust):  # pylint: disable=unused-argument
+def test_initial_solution(
+    mock_show, example_plain_env, calisto_robust
+):  # pylint: disable=unused-argument
     """Tests the initial_solution option of the Flight class. This test simply
     simulates the flight using the initial_solution option and checks if the
     all_info method returns None.
@@ -339,7 +351,9 @@ def test_initial_solution(mock_show, example_plain_env, calisto_robust):  # pyli
 
 
 @patch("matplotlib.pyplot.show")
-def test_empty_motor_flight(mock_show, example_plain_env, calisto_motorless):  # pylint: disable=unused-argument
+def test_empty_motor_flight(
+    mock_show, example_plain_env, calisto_motorless
+):  # pylint: disable=unused-argument
     flight = Flight(
         rocket=calisto_motorless,
         environment=example_plain_env,
@@ -438,3 +452,41 @@ def test_rocket_csys_equivalence(
         flight_calisto_robust.initial_solution,
         flight_calisto_nose_to_tail_robust.initial_solution,
     )
+
+
+def test_opening_shock_recorded_during_flight(calisto, example_plain_env):
+    """
+    Testing if the opening shock is being saved correctly during simulations.
+    """
+    # Defining test parachute
+    calisto.parachutes = []
+
+    target_coeff = 1.75
+    main_chute = Parachute(
+        name="Main Test",
+        cd_s=5.0,
+        trigger="apogee",
+        sampling_rate=100,
+        opening_shock_coefficient=target_coeff,
+    )
+
+    calisto.parachutes.append(main_chute)
+
+    # Simulating
+    flight = Flight(
+        rocket=calisto,
+        environment=example_plain_env,
+        rail_length=5,
+        inclination=85,
+        heading=0,
+        terminate_on_apogee=False,
+    )
+
+    # Analysing results
+    assert len(flight.parachute_events) > 0, "No parachute event registered!"
+
+    event_time, flown_chute = flight.parachute_events[0]
+
+    assert flown_chute.opening_shock_force is not None
+    assert flown_chute.opening_shock_force > 0
+    assert flown_chute.opening_shock_coefficient == target_coeff

tests/unit/test_parachute_shock.py

@@ -0,0 +1,40 @@
+from rocketpy import Parachute
+import pytest
+
+
+# TODO: Analyse if, instead of a new file, test could be included in a existing one
+def test_knacke_opening_shock_example():
+    """
+    Verifies the opening shock force calculation comparing to a textbook example.
+
+    Reference: Knacke, T. W. (1992). Parachute Recovery Systems Design Manual.
+    (Page 5-51, Figure 5-21)
+    """
+    # Setup
+    knacke_cd = 0.49
+    knacke_s = 17.76  # m^2
+    knacke_cx = 1.088
+    knacke_x1 = 1.0
+
+    knacke_density = 0.458  # kg/m^3
+    knacke_velocity = 123.2  # m/s
+
+    # Expected result
+    knacke_force = 32916.8  # N
+
+    # Defining example parachute
+    parachute = Parachute(
+        name="B-47 Test Chute",
+        cd_s=knacke_cd * knacke_s,
+        trigger="apogee",
+        sampling_rate=100,
+        opening_shock_coefficient=knacke_cx * knacke_x1,
+    )
+
+    # Calculating the shock force
+    calculated_force = parachute.calculate_opening_shock(
+        knacke_density, knacke_velocity
+    )
+
+    # Analysing results
+    assert calculated_force == pytest.approx(knacke_force, rel=1e-2)