diff --git a/CHANGELOG.md b/CHANGELOG.md index 074e28af0..507bf40b9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,6 +31,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). Attention: The newest changes should be on top --> ### Added +- ENH: Enable only radial burning [#815](https://github.com/RocketPy-Team/RocketPy/pull/815) - ENH: Tank Fluids with Variable Density from Temperature and Pressure [#852](https://github.com/RocketPy-Team/RocketPy/pull/852) - ENH: Controller (AirBrakes) and Sensors Encoding [#849](https://github.com/RocketPy-Team/RocketPy/pull/849) - EHN: Addition of ensemble variable to ECMWF dictionaries [#842](https://github.com/RocketPy-Team/RocketPy/pull/842) @@ -476,4 +477,4 @@ You can install this version by running `pip install rocketpy==1.0.1` ### Fixed - BUG: Remove NoseCone Warning [#428](https://github.com/RocketPy-Team/RocketPy/pull/428) -- BUG: motor coordinates [#423](https://github.com/RocketPy-Team/RocketPy/pull/423) +- BUG: motor coordinates [#423](https://github.com/RocketPy-Team/RocketPy/pull/423) \ No newline at end of file diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index 7cb28670c..c72af6488 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -193,8 +193,12 @@ class HybridMotor(Motor): HybridMotor.reference_pressure : int, float Atmospheric pressure in Pa at which the thrust data was recorded. It will allow to obtain the net thrust in the Flight class. + SolidMotor.only_radial_burn : bool + If True, grain regression is restricted to radial burn only (inner radius growth). + Grain length remains constant throughout the burn. Default is False. """ + # pylint: disable=too-many-arguments def __init__( # pylint: disable=too-many-arguments self, thrust_source, @@ -216,6 +220,7 @@ def __init__( # pylint: disable=too-many-arguments interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, + only_radial_burn=True, ): """Initialize Motor class, process thrust curve and geometrical parameters and store results. @@ -313,6 +318,11 @@ class Function. Thrust units are Newtons. "nozzle_to_combustion_chamber". reference_pressure : int, float, optional Atmospheric pressure in Pa at which the thrust data was recorded. + only_radial_burn : boolean, optional + If True, inhibits the grain from burning axially, only computing + radial burn. If False, allows the grain to also burn + axially. May be useful for axially inhibited grains or hybrid motors. + Default is False. Returns ------- @@ -364,6 +374,7 @@ class Function. Thrust units are Newtons. interpolation_method, coordinate_system_orientation, reference_pressure, + only_radial_burn, ) self.positioned_tanks = self.liquid.positioned_tanks diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py index 8a00eeec9..ae5d95c76 100644 --- a/rocketpy/motors/solid_motor.py +++ b/rocketpy/motors/solid_motor.py @@ -193,6 +193,9 @@ class SolidMotor(Motor): SolidMotor.reference_pressure : int, float Atmospheric pressure in Pa at which the thrust data was recorded. It will allow to obtain the net thrust in the Flight class. + SolidMotor.only_radial_burn : bool + If True, grain regression is restricted to radial burn only (inner radius growth). + Grain length remains constant throughout the burn. Default is False. """ # pylint: disable=too-many-arguments @@ -217,6 +220,7 @@ def __init__( interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, + only_radial_burn=False, ): """Initialize Motor class, process thrust curve and geometrical parameters and store results. @@ -314,6 +318,11 @@ class Function. Thrust units are Newtons. "nozzle_to_combustion_chamber". reference_pressure : int, float, optional Atmospheric pressure in Pa at which the thrust data was recorded. + only_radial_burn : boolean, optional + If True, inhibits the grain from burning axially, only computing + radial burn. If False, allows the grain to also burn + axially. May be useful for axially inhibited grains or hybrid motors. + Default is False. Returns ------- @@ -353,6 +362,9 @@ class Function. Thrust units are Newtons. ) self.grain_initial_mass = self.grain_density * self.grain_initial_volume + # Burn surface definition + self.only_radial_burn = only_radial_burn + self.evaluate_geometry() # Initialize plots and prints object @@ -405,6 +417,10 @@ def exhaust_velocity(self): self.total_impulse / self.propellant_initial_mass ).set_discrete_based_on_model(self.thrust) + # @property + # def only_radial_burn(self): + # return self._only_radial_burn + @property def propellant_initial_mass(self): """Returns the initial propellant mass. @@ -500,17 +516,25 @@ def geometry_dot(t, y): # Compute state vector derivative grain_inner_radius, grain_height = y - burn_area = ( - 2 - * np.pi - * ( - grain_outer_radius**2 - - grain_inner_radius**2 - + grain_inner_radius * grain_height + if self.only_radial_burn: + burn_area = 2 * np.pi * (grain_inner_radius * grain_height) + + grain_inner_radius_derivative = -volume_diff / burn_area + grain_height_derivative = 0 # Set to zero to disable axial burning + + else: + burn_area = ( + 2 + * np.pi + * ( + grain_outer_radius**2 + - grain_inner_radius**2 + + grain_inner_radius * grain_height + ) ) - ) - grain_inner_radius_derivative = -volume_diff / burn_area - grain_height_derivative = -2 * grain_inner_radius_derivative + + grain_inner_radius_derivative = -volume_diff / burn_area + grain_height_derivative = -2 * grain_inner_radius_derivative return [grain_inner_radius_derivative, grain_height_derivative] @@ -521,32 +545,55 @@ def geometry_jacobian(t, y): # Compute jacobian grain_inner_radius, grain_height = y - factor = volume_diff / ( - 2 - * np.pi - * ( - grain_outer_radius**2 - - grain_inner_radius**2 - + grain_inner_radius * grain_height + if self.only_radial_burn: + factor = volume_diff / ( + 2 * np.pi * (grain_inner_radius * grain_height) ** 2 ) - ** 2 - ) - inner_radius_derivative_wrt_inner_radius = factor * ( - grain_height - 2 * grain_inner_radius - ) - inner_radius_derivative_wrt_height = factor * grain_inner_radius - height_derivative_wrt_inner_radius = ( - -2 * inner_radius_derivative_wrt_inner_radius - ) - height_derivative_wrt_height = -2 * inner_radius_derivative_wrt_height - return [ - [ - inner_radius_derivative_wrt_inner_radius, - inner_radius_derivative_wrt_height, - ], - [height_derivative_wrt_inner_radius, height_derivative_wrt_height], - ] + inner_radius_derivative_wrt_inner_radius = factor * ( + grain_height - 2 * grain_inner_radius + ) + inner_radius_derivative_wrt_height = 0 + height_derivative_wrt_inner_radius = 0 + height_derivative_wrt_height = 0 + # Height is a constant, so all the derivatives with respect to it are set to zero + + return [ + [ + inner_radius_derivative_wrt_inner_radius, + inner_radius_derivative_wrt_height, + ], + [height_derivative_wrt_inner_radius, height_derivative_wrt_height], + ] + + else: + factor = volume_diff / ( + 2 + * np.pi + * ( + grain_outer_radius**2 + - grain_inner_radius**2 + + grain_inner_radius * grain_height + ) + ** 2 + ) + + inner_radius_derivative_wrt_inner_radius = factor * ( + grain_height - 2 * grain_inner_radius + ) + inner_radius_derivative_wrt_height = factor * grain_inner_radius + height_derivative_wrt_inner_radius = ( + -2 * inner_radius_derivative_wrt_inner_radius + ) + height_derivative_wrt_height = -2 * inner_radius_derivative_wrt_height + + return [ + [ + inner_radius_derivative_wrt_inner_radius, + inner_radius_derivative_wrt_height, + ], + [height_derivative_wrt_inner_radius, height_derivative_wrt_height], + ] def terminate_burn(t, y): # pylint: disable=unused-argument end_function = (self.grain_outer_radius - y[0]) * y[1] @@ -597,16 +644,24 @@ def burn_area(self): burn_area : Function Function representing the burn area progression with the time. """ - burn_area = ( - 2 - * np.pi - * ( - self.grain_outer_radius**2 - - self.grain_inner_radius**2 - + self.grain_inner_radius * self.grain_height + if self.only_radial_burn: + burn_area = ( + 2 + * np.pi + * (self.grain_inner_radius * self.grain_height) + * self.grain_number + ) + else: + burn_area = ( + 2 + * np.pi + * ( + self.grain_outer_radius**2 + - self.grain_inner_radius**2 + + self.grain_inner_radius * self.grain_height + ) + * self.grain_number ) - * self.grain_number - ) return burn_area @funcify_method("Time (s)", "burn rate (m/s)") @@ -778,6 +833,7 @@ def to_dict(self, **kwargs): "grain_initial_height": self.grain_initial_height, "grain_separation": self.grain_separation, "grains_center_of_mass_position": self.grains_center_of_mass_position, + "only_radial_burn": self.only_radial_burn, } ) @@ -827,4 +883,5 @@ def from_dict(cls, data): interpolation_method=data["interpolate"], coordinate_system_orientation=data["coordinate_system_orientation"], reference_pressure=data.get("reference_pressure"), + only_radial_burn=data.get("only_radial_burn", False), ) diff --git a/tests/fixtures/motor/hybrid_fixtures.py b/tests/fixtures/motor/hybrid_fixtures.py index 923a640b1..35812fbfb 100644 --- a/tests/fixtures/motor/hybrid_fixtures.py +++ b/tests/fixtures/motor/hybrid_fixtures.py @@ -4,7 +4,7 @@ @pytest.fixture -def hybrid_motor(spherical_oxidizer_tank): +def hybrid_motor(oxidizer_tank): """An example of a hybrid motor with spherical oxidizer tank and fuel grains. @@ -35,6 +35,6 @@ def hybrid_motor(spherical_oxidizer_tank): grains_center_of_mass_position=-0.1, ) - motor.add_tank(spherical_oxidizer_tank, position=0.3) + motor.add_tank(oxidizer_tank, position=0.3) return motor diff --git a/tests/integration/test_solidmotor.py b/tests/integration/test_solidmotor.py new file mode 100644 index 000000000..a9aef6a8b --- /dev/null +++ b/tests/integration/test_solidmotor.py @@ -0,0 +1,16 @@ +from unittest.mock import patch + + +@patch("matplotlib.pyplot.show") +def test_solid_motor_info(mock_show, cesaroni_m1670): # pylint: disable=unused-argument + """Tests the SolidMotor.all_info() method. + + Parameters + ---------- + mock_show : mock + Mock of the matplotlib.pyplot.show function. + cesaroni_m1670 : rocketpy.SolidMotor + The SolidMotor object to be used in the tests. + """ + assert cesaroni_m1670.info() is None + assert cesaroni_m1670.all_info() is None diff --git a/tests/unit/test_hybridmotor.py b/tests/unit/test_hybridmotor.py index ef03a1998..e6b69d555 100644 --- a/tests/unit/test_hybridmotor.py +++ b/tests/unit/test_hybridmotor.py @@ -56,7 +56,7 @@ def test_hybrid_motor_basic_parameters(hybrid_motor): assert hybrid_motor.liquid.positioned_tanks[0]["position"] == 0.3 -def test_hybrid_motor_thrust_parameters(hybrid_motor, spherical_oxidizer_tank): +def test_hybrid_motor_thrust_parameters(hybrid_motor, oxidizer_tank): """Tests the HybridMotor class thrust parameters. Parameters @@ -77,13 +77,13 @@ def test_hybrid_motor_thrust_parameters(hybrid_motor, spherical_oxidizer_tank): * GRAIN_INITIAL_HEIGHT * GRAIN_NUMBER ) - initial_oxidizer_mass = spherical_oxidizer_tank.fluid_mass(0) + initial_oxidizer_mass = oxidizer_tank.fluid_mass(0) initial_mass = initial_grain_mass + initial_oxidizer_mass expected_exhaust_velocity = expected_total_impulse / initial_mass expected_mass_flow_rate = -expected_thrust / expected_exhaust_velocity expected_grain_mass_flow_rate = ( - expected_mass_flow_rate - spherical_oxidizer_tank.net_mass_flow_rate + expected_mass_flow_rate - oxidizer_tank.net_mass_flow_rate ) assert pytest.approx(hybrid_motor.thrust.y_array) == expected_thrust.y_array @@ -100,7 +100,7 @@ def test_hybrid_motor_thrust_parameters(hybrid_motor, spherical_oxidizer_tank): ) == expected_grain_mass_flow_rate(t) -def test_hybrid_motor_center_of_mass(hybrid_motor, spherical_oxidizer_tank): +def test_hybrid_motor_center_of_mass(hybrid_motor, oxidizer_tank): """Tests the HybridMotor class center of mass. Parameters @@ -110,25 +110,25 @@ def test_hybrid_motor_center_of_mass(hybrid_motor, spherical_oxidizer_tank): spherical_oxidizer_tank : rocketpy.SphericalTank The SphericalTank object to be used in the tests. """ - oxidizer_mass = spherical_oxidizer_tank.fluid_mass + oxidizer_mass = oxidizer_tank.fluid_mass grain_mass = hybrid_motor.solid.propellant_mass propellant_balance = grain_mass * GRAINS_CENTER_OF_MASS_POSITION + oxidizer_mass * ( - OXIDIZER_TANK_POSITION + spherical_oxidizer_tank.center_of_mass + OXIDIZER_TANK_POSITION + oxidizer_tank.center_of_mass ) balance = propellant_balance + DRY_MASS * CENTER_OF_DRY_MASS propellant_center_of_mass = propellant_balance / (grain_mass + oxidizer_mass) center_of_mass = balance / (grain_mass + oxidizer_mass + DRY_MASS) - for t in np.linspace(0, 100, 100): + for t in np.linspace(0, BURN_TIME, 100): assert pytest.approx( hybrid_motor.center_of_propellant_mass(t) ) == propellant_center_of_mass(t) assert pytest.approx(hybrid_motor.center_of_mass(t)) == center_of_mass(t) -def test_hybrid_motor_inertia(hybrid_motor, spherical_oxidizer_tank): +def test_hybrid_motor_inertia(hybrid_motor, oxidizer_tank): """Tests the HybridMotor class inertia. Parameters @@ -138,8 +138,8 @@ def test_hybrid_motor_inertia(hybrid_motor, spherical_oxidizer_tank): spherical_oxidizer_tank : rocketpy.SphericalTank The SphericalTank object to be used in the tests. """ - oxidizer_mass = spherical_oxidizer_tank.fluid_mass - oxidizer_inertia = spherical_oxidizer_tank.inertia + oxidizer_mass = oxidizer_tank.fluid_mass + oxidizer_inertia = oxidizer_tank.inertia grain_mass = hybrid_motor.solid.propellant_mass grain_inertia = hybrid_motor.solid.propellant_I_11 propellant_mass = oxidizer_mass + grain_mass @@ -153,7 +153,7 @@ def test_hybrid_motor_inertia(hybrid_motor, spherical_oxidizer_tank): oxidizer_mass * ( OXIDIZER_TANK_POSITION - + spherical_oxidizer_tank.center_of_mass + + oxidizer_tank.center_of_mass - hybrid_motor.center_of_propellant_mass ) ** 2 @@ -170,9 +170,46 @@ def test_hybrid_motor_inertia(hybrid_motor, spherical_oxidizer_tank): + DRY_MASS * (-hybrid_motor.center_of_mass + CENTER_OF_DRY_MASS) ** 2 ) - for t in np.linspace(0, 100, 100): + for t in np.linspace(0, BURN_TIME, 100): assert pytest.approx(hybrid_motor.propellant_I_11(t)) == propellant_inertia(t) assert pytest.approx(hybrid_motor.I_11(t)) == inertia(t) # Assert cylindrical symmetry assert pytest.approx(hybrid_motor.propellant_I_22(t)) == propellant_inertia(t) + + +def test_hybrid_motor_only_radial_burn_behavior(hybrid_motor): + """ + Test if only_radial_burn flag in HybridMotor propagates to its SolidMotor + and affects burn_area calculation. + """ + motor = hybrid_motor + + # Activates the radial burning + motor.solid.only_radial_burn = True + assert motor.solid.only_radial_burn is True + + # Calculates the expected initial area + burn_area_radial = ( + 2 + * np.pi + * (motor.solid.grain_inner_radius(0) * motor.solid.grain_height(0)) + * motor.solid.grain_number + ) + + assert np.isclose(motor.solid.burn_area(0), burn_area_radial, atol=1e-12) + + # Deactivates the radial burning and recalculate the geometry + motor.solid.only_radial_burn = False + motor.solid.evaluate_geometry() + assert motor.solid.only_radial_burn is False + + # In this case the burning area also considers the bases of the grain + inner_radius = motor.solid.grain_inner_radius(0) + outer_radius = motor.solid.grain_outer_radius + burn_area_total = ( + burn_area_radial + + 2 * np.pi * (outer_radius**2 - inner_radius**2) * motor.solid.grain_number + ) + assert np.isclose(motor.solid.burn_area(0), burn_area_total, atol=1e-12) + assert motor.solid.burn_area(0) > burn_area_radial diff --git a/tests/unit/test_solidmotor.py b/tests/unit/test_solidmotor.py index 3f829d222..4f61a2dc2 100644 --- a/tests/unit/test_solidmotor.py +++ b/tests/unit/test_solidmotor.py @@ -279,3 +279,59 @@ def test_reshape_thrust_curve_asserts_resultant_thrust_curve_correct( assert thrust_reshaped[1][1] == 100 * (tuple_parametric[1] / 7539.1875) assert thrust_reshaped[7][1] == 2034 * (tuple_parametric[1] / 7539.1875) + + +def test_only_radial_burn_parameter_effect(cesaroni_m1670): + """Tests the effect of the only_radial_burn parameter on burn area + calculation. When enabled, the burn area should only account for + the radial surface of the grains (no axial regression). + + Parameters + ---------- + cesaroni_m1670 : rocketpy.SolidMotor + The SolidMotor object used in the test. + """ + motor = cesaroni_m1670 + motor.only_radial_burn = True + assert motor.only_radial_burn is True + + # When only_radial_burn is active, burn_area should consider only radial area + burn_area_radial = ( + 2 + * np.pi + * (motor.grain_inner_radius(0) * motor.grain_height(0)) + * motor.grain_number + ) + assert np.isclose(motor.burn_area(0), burn_area_radial, atol=1e-12) + + +def test_evaluate_geometry_updates_properties(cesaroni_m1670): + """Tests if the grain geometry evaluation correctly updates SolidMotor + properties after instantiation. It ensures that grain geometry + functions are created and behave as expected. + + Parameters + ---------- + cesaroni_m1670 : rocketpy.SolidMotor + The SolidMotor object used in the test. + """ + motor = cesaroni_m1670 + + assert hasattr(motor, "grain_inner_radius") + assert hasattr(motor, "grain_height") + + # Checks if the domain of grain_inner_radius function is consistent + times = motor.grain_inner_radius.x_array + values = motor.grain_inner_radius.y_array + + assert times[0] == 0 # expected initial time + assert ( + values[0] == motor.grain_initial_inner_radius + ) # expected initial inner radius + assert ( + values[-1] <= motor.grain_outer_radius + ) # final inner radius should be less or equal than outer radius + + # Evaluates without error + val = motor.grain_inner_radius(0.5) # evaluate at intermediate time + assert isinstance(val, float)