From 531cf381d2c05666a167023fe2bf938e1da5757d Mon Sep 17 00:00:00 2001 From: Felipe Almeida Date: Sat, 6 Dec 2025 11:50:12 -0300 Subject: [PATCH 1/2] ENH: implement parachute opening shock force estimation --- rocketpy/rocket/parachute.py | 36 +++++++++++ rocketpy/simulation/flight.py | 46 +++++++++----- tests/integration/simulation/test_flight.py | 69 ++++++++++++++++++--- tests/unit/test_parachute_shock.py | 40 ++++++++++++ 4 files changed, 166 insertions(+), 25 deletions(-) create mode 100644 tests/unit/test_parachute_shock.py diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index f0bf86f67..82739be6c 100644 --- a/rocketpy/rocket/parachute.py +++ b/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 @@ -346,3 +357,28 @@ def from_dict(cls, data): ) 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 diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 3e0204ece..27c91143e 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -773,6 +773,22 @@ def __simulate(self, verbose): self.y_sol, self.sensors, ): + + # 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 + ) + ) + # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create phase for time after detection and before inflation @@ -800,8 +816,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 +1063,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, @@ -1601,7 +1611,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 +2159,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 +3996,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) diff --git a/tests/integration/simulation/test_flight.py b/tests/integration/simulation/test_flight.py index f40eb6b27..3be176a9f 100644 --- a/tests/integration/simulation/test_flight.py +++ b/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,42 @@ def test_rocket_csys_equivalence( flight_calisto_robust.initial_solution, flight_calisto_nose_to_tail_robust.initial_solution, ) + + +# TODO: fix the issues on this test and debug shock analysis +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 diff --git a/tests/unit/test_parachute_shock.py b/tests/unit/test_parachute_shock.py new file mode 100644 index 000000000..bda3f2322 --- /dev/null +++ b/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) From 0f8f3a66df0e87eb29aac5ec02db401cba6bdb16 Mon Sep 17 00:00:00 2001 From: Felipe Almeida Date: Sun, 7 Dec 2025 11:12:08 -0300 Subject: [PATCH 2/2] Fix: Add opening shock calculation to overshoot loop --- rocketpy/rocket/parachute.py | 3 ++ rocketpy/simulation/flight.py | 37 +++++++++++++-------- tests/integration/simulation/test_flight.py | 1 - 3 files changed, 26 insertions(+), 15 deletions(-) diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index 82739be6c..98b651b67 100644 --- a/rocketpy/rocket/parachute.py +++ b/rocketpy/rocket/parachute.py @@ -321,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): @@ -332,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 @@ -354,6 +356,7 @@ 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 diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 27c91143e..c9d1c13fd 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -774,20 +774,8 @@ def __simulate(self, verbose): self.sensors, ): - # 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 - ) - ) + # Calculates the parachute's opening shock force + self.calculate_parachute_opening_shock_force(parachute) # Remove parachute from flight parachutes self.parachutes.remove(parachute) @@ -1108,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] @@ -4371,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 + ) diff --git a/tests/integration/simulation/test_flight.py b/tests/integration/simulation/test_flight.py index 3be176a9f..093a57d8e 100644 --- a/tests/integration/simulation/test_flight.py +++ b/tests/integration/simulation/test_flight.py @@ -454,7 +454,6 @@ def test_rocket_csys_equivalence( ) -# TODO: fix the issues on this test and debug shock analysis def test_opening_shock_recorded_during_flight(calisto, example_plain_env): """ Testing if the opening shock is being saved correctly during simulations.