diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index f0bf86f67..98b651b67 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 @@ -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 diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 3e0204ece..c9d1c13fd 100644 --- a/rocketpy/simulation/flight.py +++ b/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 + ) diff --git a/tests/integration/simulation/test_flight.py b/tests/integration/simulation/test_flight.py index f40eb6b27..093a57d8e 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,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 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)