ENH: estimate parachute opening shock force (#161)

Author: wuisabel-gif

rocketpy/rocket/parachute.py

@@ -123,6 +123,11 @@ class Parachute:
     Parachute.added_mass_coefficient : float
         Coefficient used to calculate the added-mass due to dragged air. It is
         calculated from the porosity of the parachute.
+    Parachute.opening_shock_coefficient : float
+        Dimensionless opening-force coefficient (``Cx``) used to estimate the
+        peak opening shock force via
+        :meth:`Parachute.evaluate_opening_shock_force`. ``None`` when the
+        estimate is disabled.
     """
 
     def __init__(
@@ -137,6 +142,7 @@ def __init__(
         height=None,
         porosity=0.0432,
         drag_coefficient=1.4,
+        opening_shock_coefficient=None,
     ):
         """Initializes Parachute class.
 
@@ -217,6 +223,15 @@ def __init__(
             - **1.5** — extended-skirt canopy
 
             Has no effect when ``radius`` is explicitly provided.
+        opening_shock_coefficient : float, optional
+            Dimensionless opening-force coefficient (``Cx``) used to estimate
+            the peak parachute opening shock force via
+            :meth:`evaluate_opening_shock_force`, following Knacke's
+            *Parachute Recovery Systems Design Manual*. It lumps together the
+            opening-force coefficient and the infinite-mass factor (``X1``).
+            Typical values range from about 1.2 to 2.0 depending on canopy
+            type. If ``None`` (default), the opening shock force is not
+            estimated. Units are dimensionless.
         """
 
         # Save arguments as attributes
@@ -228,6 +243,7 @@ def __init__(
         self.noise = noise
         self.drag_coefficient = drag_coefficient
         self.porosity = porosity
+        self.opening_shock_coefficient = opening_shock_coefficient
 
         # Initialize derived attributes
         self.radius = self.__resolve_radius(radius, cd_s, drag_coefficient)
@@ -259,6 +275,43 @@ def __compute_added_mass_coefficient(self, porosity):
             1 - 1.465 * porosity - 0.25975 * porosity**2 + 1.2626 * porosity**3
         )
 
+    def evaluate_opening_shock_force(self, air_density, velocity):
+        """Estimate the peak parachute opening shock force.
+
+        Uses the standard approximation from Knacke's *Parachute Recovery
+        Systems Design Manual*: ``F_o = Cx * q * cd_s``, where ``q`` is the
+        dynamic pressure ``0.5 * air_density * velocity ** 2`` at inflation and
+        ``Cx`` is the ``opening_shock_coefficient``. This is an empirical
+        estimate of the peak transient (inflation) load on the recovery
+        harness, not a result of the equation-of-motion integration, and it is
+        typically several times larger than the steady-state drag force
+        ``q * cd_s``.
+
+        Parameters
+        ----------
+        air_density : float
+            Air density at inflation, in kg/m^3.
+        velocity : float
+            Freestream speed at inflation, in m/s.
+
+        Returns
+        -------
+        float
+            Peak opening shock force, in newtons.
+
+        Raises
+        ------
+        ValueError
+            If ``opening_shock_coefficient`` was not set at construction.
+        """
+        if self.opening_shock_coefficient is None:
+            raise ValueError(
+                "opening_shock_coefficient must be set to estimate the "
+                "opening shock force."
+            )
+        dynamic_pressure = 0.5 * air_density * velocity**2
+        return self.opening_shock_coefficient * dynamic_pressure * self.cd_s
+
     def __init_noise(self, noise):
         """Initializes all noise-related attributes.
 
@@ -421,6 +474,7 @@ def to_dict(self, **kwargs):
             "drag_coefficient": self.drag_coefficient,
             "height": self.height,
             "porosity": self.porosity,
+            "opening_shock_coefficient": self.opening_shock_coefficient,
         }
 
         if kwargs.get("include_outputs", False):
@@ -455,6 +509,7 @@ def from_dict(cls, data):
             drag_coefficient=data.get("drag_coefficient", 1.4),
             height=data.get("height", None),
             porosity=data.get("porosity", 0.0432),
+            opening_shock_coefficient=data.get("opening_shock_coefficient", None),
         )
 
         return parachute

tests/unit/rocket/test_parachute.py

@@ -109,3 +109,58 @@ def test_from_dict_defaults_drag_coefficient_to_1_4_when_absent(self):
         }
         parachute = Parachute.from_dict(data)
         assert parachute.drag_coefficient == pytest.approx(1.4)
+
+
+class TestParachuteOpeningShockForce:
+    """Tests for the opening shock force estimation (issue #161)."""
+
+    def test_opening_shock_force_matches_knacke_formula(self):
+        """The estimate must equal Cx * q * cd_s with q = 0.5 * rho * v**2."""
+        cd_s = 10.0
+        cx = 1.6
+        air_density = 1.05
+        velocity = 30.0
+        parachute = _make_parachute(cd_s=cd_s, opening_shock_coefficient=cx)
+
+        force = parachute.evaluate_opening_shock_force(air_density, velocity)
+
+        expected = cx * (0.5 * air_density * velocity**2) * cd_s
+        assert force == pytest.approx(expected, rel=1e-12)
+
+    def test_opening_shock_force_exceeds_steady_drag(self):
+        """For a coefficient above 1, the peak load must exceed the steady
+        drag force q * cd_s at the same condition."""
+        cd_s = 10.0
+        air_density = 1.05
+        velocity = 30.0
+        parachute = _make_parachute(cd_s=cd_s, opening_shock_coefficient=1.6)
+
+        force = parachute.evaluate_opening_shock_force(air_density, velocity)
+
+        steady_drag = (0.5 * air_density * velocity**2) * cd_s
+        assert force > steady_drag
+
+    def test_opening_shock_force_scales_with_velocity_squared(self):
+        """Doubling the inflation velocity must quadruple the peak load."""
+        parachute = _make_parachute(opening_shock_coefficient=1.5)
+
+        force_low = parachute.evaluate_opening_shock_force(1.2, 20.0)
+        force_high = parachute.evaluate_opening_shock_force(1.2, 40.0)
+
+        assert force_high == pytest.approx(4.0 * force_low, rel=1e-12)
+
+    def test_opening_shock_force_requires_coefficient(self):
+        """Without an opening_shock_coefficient the estimate is disabled."""
+        parachute = _make_parachute()  # opening_shock_coefficient defaults to None
+
+        assert parachute.opening_shock_coefficient is None
+        with pytest.raises(ValueError, match="opening_shock_coefficient"):
+            parachute.evaluate_opening_shock_force(1.2, 30.0)
+
+    def test_opening_shock_coefficient_survives_dict_round_trip(self):
+        """to_dict/from_dict must preserve the opening_shock_coefficient."""
+        parachute = _make_parachute(opening_shock_coefficient=1.7)
+
+        restored = Parachute.from_dict(parachute.to_dict())
+
+        assert restored.opening_shock_coefficient == 1.7