ENH: Parachute Inflation Modeling

rocketpy/rocket/parachute.py

@@ -1,5 +1,6 @@
 from inspect import signature
 
+import math
 import numpy as np
 
 from rocketpy.tools import from_hex_decode, to_hex_encode
@@ -98,6 +99,9 @@ class Parachute:
     Parachute.height : float, None
         Length of the unique semi-axis (height) of the inflated hemispheroid
         parachute in meters.
+    Parachute.initial_radius : float, None
+        Initial radius of the parachute on deployment in meters. Used to model the
+        parachute inflation.
     Parachute.porosity : float
         Geometric porosity of the canopy (ratio of open area to total canopy area),
         in [0, 1]. Affects only the added-mass scaling during descent; it does
@@ -118,6 +122,7 @@ def __init__(
         noise=(0, 0, 0),
         radius=1.5,
         height=None,
+        initial_radius=None,
         porosity=0.0432,
     ):
         """Initializes Parachute class.
@@ -179,6 +184,9 @@ def __init__(
             Length of the unique semi-axis (height) of the inflated hemispheroid
             parachute. Default value is the radius of the parachute.
             Units are in meters.
+        initial_radius : float, optional
+            Initial radius of the parachute on deployment in meters. Used to model the
+            parachute inflation. Default value is the parachute radius (no inflation).
         porosity : float, optional
             Geometric porosity of the canopy (ratio of open area to total canopy area),
             in [0, 1]. Affects only the added-mass scaling during descent; it does
@@ -202,14 +210,20 @@ def __init__(
         self.noise_signal_function = Function(0)
         self.radius = radius
         self.height = height or radius
+        self.initial_radius = initial_radius or radius
+        self.initial_volume = (
+            (4 / 3)
+            * math.pi
+            * (self.parachute_height / self.parachute_radius)
+            * (min(self.parachute_radius, self.rocket.radius)) ** 3
+        )
         self.porosity = porosity
         self.added_mass_coefficient = 1.068 * (
             1
             - 1.465 * self.porosity
             - 0.25975 * self.porosity**2
             + 1.2626 * self.porosity**3
         )
-
         alpha, beta = self.noise_corr
         self.noise_function = lambda: alpha * self.noise_signal[-1][
             1
@@ -309,6 +323,7 @@ def to_dict(self, **kwargs):
             "noise": self.noise,
             "radius": self.radius,
             "height": self.height,
+            "initial_radius": self.initial_radius,
             "porosity": self.porosity,
         }
 
@@ -342,6 +357,7 @@ def from_dict(cls, data):
             noise=data["noise"],
             radius=data.get("radius", 1.5),
             height=data.get("height", None),
+            initial_radius=data.get("initial_radius", None),
             porosity=data.get("porosity", 0.0432),
         )
 

rocketpy/rocket/rocket.py

@@ -1485,6 +1485,7 @@ def add_parachute(
         noise=(0, 0, 0),
         radius=1.5,
         height=None,
+        initial_radius=None,
         porosity=0.0432,
     ):
         """Creates a new parachute, storing its parameters such as
@@ -1552,6 +1553,9 @@ def add_parachute(
             Length of the unique semi-axis (height) of the inflated hemispheroid
             parachute. Default value is the radius of the parachute.
             Units are in meters.
+        initial_radius : float, optional
+            Initial radius of the parachute on deployment in meters. Used to model the
+            parachute inflation. Default value is the parachute radius (no inflation).
         porosity : float, optional
             Geometric porosity of the canopy (ratio of open area to total canopy area),
             in [0, 1]. Affects only the added-mass scaling during descent; it does
@@ -1567,15 +1571,16 @@ def add_parachute(
             Flight simulation.
         """
         parachute = Parachute(
-            name,
-            cd_s,
-            trigger,
-            sampling_rate,
-            lag,
-            noise,
-            radius,
-            height,
-            porosity,
+            name=name,
+            cd_s=cd_s,
+            trigger=trigger,
+            sampling_rate=sampling_rate,
+            lag=lag,
+            noise=noise,
+            radius=radius,
+            height=height,
+            initial_radius=initial_radius,
+            porosity=porosity,
         )
         self.parachutes.append(parachute)
         return self.parachutes[-1]

rocketpy/simulation/flight.py

@@ -765,6 +765,15 @@ def __simulate(self, verbose):
                                 "parachute_added_mass_coefficient",
                                 added_mass_coefficient,
                             ),
+                            lambda self,
+                            initial_volume=parachute.initial_volume: setattr(
+                                self,
+                                "parachute_volume",
+                                initial_volume,
+                            ),
+                            lambda self: delattr(self, "__t0")
+                            if hasattr(self, "__t0")
+                            else None,
                         ]
                         self.flight_phases.add_phase(
                             node.t + parachute.lag,
@@ -987,6 +996,12 @@ def __check_and_handle_parachute_triggers(
                     "parachute_added_mass_coefficient",
                     added_mass_coefficient,
                 ),
+                lambda self, initial_volume=parachute.initial_volume: setattr(
+                    self,
+                    "parachute_volume",
+                    initial_volume,
+                ),
+                lambda self: delattr(self, "__t0") if hasattr(self, "__t0") else None,
             ]
             self.flight_phases.add_phase(
                 node.t + parachute.lag,
@@ -2702,31 +2717,58 @@ def u_dot_parachute(self, t, u, post_processing=False):
         # Get the mass of the rocket
         mp = self.rocket.dry_mass
 
-        # to = 1.2
-        # eta = 1
-        # Rdot = (6 * R * (1 - eta) / (1.2**6)) * (
-        #     (1 - eta) * t**5 + eta * (to**3) * (t**2)
-        # )
-        # Rdot = 0
-
-        # tf = 8 * nominal diameter / velocity at line stretch
-
-        # Calculate added mass
-        ma = (
-            self.parachute_added_mass_coefficient
-            * rho
-            * (2 / 3)
-            * np.pi
-            * self.parachute_radius**2
-            * self.parachute_height
-        )
-
         # Calculate freestream speed
         freestream_x = vx - wind_velocity_x
         freestream_y = vy - wind_velocity_y
         freestream_z = vz
         free_stream_speed = (freestream_x**2 + freestream_y**2 + freestream_z**2) ** 0.5
 
+        # Initialize parachute geometrical parameters
+        inflated_radius = (
+            (3 * self.parachute_volume * self.parachute_radius)
+            / (4 * math.pi * self.parachute_height)
+        ) ** (1 / 3)
+        inflated_height = (
+            inflated_radius * self.parachute_height / self.parachute_radius
+        )
+
+        # Calculate the surface area of the parachute
+        if self.parachute_radius > self.parachute_height:
+            e = math.sqrt(1 - (inflated_height**2) / (inflated_radius**2))
+            surface_area = (
+                math.pi * inflated_radius**2 * (1 + (1 - e**2) / e * math.atanh(e))
+            )
+        else:
+            e = math.sqrt(1 - (inflated_radius**2) / (inflated_height**2))
+            surface_area = (
+                math.pi
+                * inflated_radius**2
+                * (1 + inflated_height / (e * inflated_radius) * math.asin(e))
+            )
+
+        # Calculate volume flow of air into parachute
+        volume_flow = (
+            freestream_z  # considering parachute as vertical
+            * (
+                (math.pi * inflated_radius**2)
+                - (self.parachute_porosity * surface_area)
+            )
+        )
+
+        # Getting time step
+        self.__t0 = getattr(self, "t0", t)
+        t1 = t
+        dt = t1 - self.__t0
+
+        # Integrating parachute volume
+        self.parachute_volume += volume_flow * dt
+
+        # Dragged air mass
+        ma = self.parachute_volume * rho
+
+        # Moving time step
+        self.__t0 = t1
+
         # Determine drag force
         pseudo_drag = -0.5 * rho * self.parachute_cd_s * free_stream_speed
         # pseudo_drag = pseudo_drag - ka * rho * 4 * np.pi * (R**2) * Rdot