ENH: Implement only_radial_burn feature for SolidMotor and HybridMotor classes

CHANGELOG.md

@@ -32,6 +32,7 @@ Attention: The newest changes should be on top -->
 
 ### Added
 
+- ENH: Enable only radial burning [#815](https://github.com/RocketPy-Team/RocketPy/pull/815)
 
 ### Changed
 

rocketpy/motors/hybrid_motor.py

@@ -216,6 +216,7 @@ def __init__(  # pylint: disable=too-many-arguments
         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.
@@ -313,6 +314,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 +370,7 @@ class Function. Thrust units are Newtons.
             interpolation_method,
             coordinate_system_orientation,
             reference_pressure,
+            only_radial_burn,
         )
 
         self.positioned_tanks = self.liquid.positioned_tanks

rocketpy/motors/solid_motor.py

@@ -217,6 +217,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 +315,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 +359,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
@@ -500,17 +509,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 +538,56 @@ 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 constant when only radial burning is enabled,
+                # so all derivatives with respect to height are 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 +638,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 +827,7 @@ def to_dict(self, include_outputs=False):
                 "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,
             }
         )
 
@@ -822,4 +872,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),
         )

tests/integration/test_environment.py

@@ -28,7 +28,9 @@ def test_set_elevation_open_elevation(
 
 
 @patch("matplotlib.pyplot.show")
-def test_era5_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_era5_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Reanalysis model with the ERA5 file. It uses an example file
     available in the data/weather folder of the RocketPy repository.
 
@@ -49,7 +51,9 @@ def test_era5_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=u
 
 
 @patch("matplotlib.pyplot.show")
-def test_custom_atmosphere(mock_show, example_plain_env):  # pylint: disable=unused-argument
+def test_custom_atmosphere(
+    mock_show, example_plain_env
+):  # pylint: disable=unused-argument
     """Tests the custom atmosphere model in the environment object.
 
     Parameters
@@ -74,7 +78,9 @@ def test_custom_atmosphere(mock_show, example_plain_env):  # pylint: disable=unu
 
 
 @patch("matplotlib.pyplot.show")
-def test_standard_atmosphere(mock_show, example_plain_env):  # pylint: disable=unused-argument
+def test_standard_atmosphere(
+    mock_show, example_plain_env
+):  # pylint: disable=unused-argument
     """Tests the standard atmosphere model in the environment object.
 
     Parameters
@@ -126,7 +132,9 @@ def test_windy_atmosphere(example_euroc_env, model_name):
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_gfs_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_gfs_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Forecast model with the GFS file. It does not test the values,
     instead the test checks if the method runs without errors.
 
@@ -143,7 +151,9 @@ def test_gfs_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=un
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_nam_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_nam_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Forecast model with the NAM file.
 
     Parameters
@@ -159,7 +169,9 @@ def test_nam_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=un
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_rap_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_rap_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     today = date.today()
     now = datetime.now(timezone.utc)
     example_spaceport_env.set_date((today.year, today.month, today.day, now.hour))
@@ -169,7 +181,9 @@ def test_rap_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=un
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_gefs_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_gefs_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Ensemble model with the GEFS file.
 
     Parameters
@@ -185,7 +199,9 @@ def test_gefs_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=u
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_wyoming_sounding_atmosphere(mock_show, example_plain_env):  # pylint: disable=unused-argument
+def test_wyoming_sounding_atmosphere(
+    mock_show, example_plain_env
+):  # pylint: disable=unused-argument
     """Asserts whether the Wyoming sounding model in the environment
     object behaves as expected with respect to some attributes such
     as pressure, barometric_height, wind_velocity and temperature.
@@ -220,7 +236,9 @@ def test_wyoming_sounding_atmosphere(mock_show, example_plain_env):  # pylint: d
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_hiresw_ensemble_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_hiresw_ensemble_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Forecast model with the HIRESW file.
 
     Parameters
@@ -246,7 +264,9 @@ def test_hiresw_ensemble_atmosphere(mock_show, example_spaceport_env):  # pylint
 
 @pytest.mark.skip(reason="CMC model is currently not working")
 @patch("matplotlib.pyplot.show")
-def test_cmc_atmosphere(mock_show, example_spaceport_env):  # pylint: disable=unused-argument
+def test_cmc_atmosphere(
+    mock_show, example_spaceport_env
+):  # pylint: disable=unused-argument
     """Tests the Ensemble model with the CMC file.
 
     Parameters

tests/integration/test_environment_analysis.py

@@ -60,5 +60,7 @@ def test_exports(mock_show, env_analysis):  # pylint: disable=unused-argument
 
 @pytest.mark.slow
 @patch("matplotlib.pyplot.show")
-def test_create_environment_object(mock_show, env_analysis):  # pylint: disable=unused-argument
+def test_create_environment_object(
+    mock_show, env_analysis
+):  # pylint: disable=unused-argument
     assert isinstance(env_analysis.create_environment_object(), Environment)