diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index f0bf86f67..0f9f039f8 100644 --- a/rocketpy/rocket/parachute.py +++ b/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,6 +210,13 @@ 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 @@ -209,7 +224,6 @@ def __init__( - 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), ) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 9ac1bb86e..5ae4f5bcc 100644 --- a/rocketpy/rocket/rocket.py +++ b/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] diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a7e0374b2..6bfbb9882 100644 --- a/rocketpy/simulation/flight.py +++ b/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