Skip to content

ENH: Enable only radial burning #815

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ Attention: The newest changes should be on top -->

### Added

- ENH: Enable only radial burning [#815](https://github.com/RocketPy-Team/RocketPy/pull/815)
- ENH: Parallel mode for monte-carlo simulations 2 [#768](https://github.com/RocketPy-Team/RocketPy/pull/768)
- DOC: ASTRA Flight Example [#770](https://github.com/RocketPy-Team/RocketPy/pull/770)
- ENH: Add Eccentricity to Stochastic Simulations [#792](https://github.com/RocketPy-Team/RocketPy/pull/792)
Expand Down
2 changes: 2 additions & 0 deletions rocketpy/motors/hybrid_motor.py
Original file line number Diff line number Diff line change
Expand Up @@ -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=True,
):
"""Initialize Motor class, process thrust curve and geometrical
parameters and store results.
Expand Down Expand Up @@ -364,6 +365,7 @@ class Function. Thrust units are Newtons.
interpolation_method,
coordinate_system_orientation,
reference_pressure,
only_radial_burn,
)

self.positioned_tanks = self.liquid.positioned_tanks
Expand Down
136 changes: 93 additions & 43 deletions rocketpy/motors/solid_motor.py
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is important to add at least an integration test to this new implementation.

Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
-------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]

Expand All @@ -521,32 +538,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]
Expand Down Expand Up @@ -597,16 +637,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)")
Expand Down Expand Up @@ -778,6 +826,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,
}
)

Expand Down Expand Up @@ -822,4 +871,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),
)
Loading