diff --git a/CHANGELOG.md b/CHANGELOG.md index 5503705ea..a91a671ec 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: Rail button bending moments calculation in Flight class [#893](https://github.com/RocketPy-Team/RocketPy/pull/893) - ENH: Built-in flight comparison tool (`FlightComparator`) to validate simulations against external data [#888](https://github.com/RocketPy-Team/RocketPy/pull/888) - ENH: Add persistent caching for ThrustCurve API [#881](https://github.com/RocketPy-Team/RocketPy/pull/881) - ENH: Compatibility with MERRA-2 atmosphere reanalysis files [#825](https://github.com/RocketPy-Team/RocketPy/pull/825) diff --git a/docs/user/flight.rst b/docs/user/flight.rst index 17dde836b..31e7ab588 100644 --- a/docs/user/flight.rst +++ b/docs/user/flight.rst @@ -266,6 +266,39 @@ The Flight object provides access to all forces and accelerations acting on the M2 = flight.M2 # Pitch moment M3 = flight.M3 # Yaw moment +Rail Button Forces and Bending Moments +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +During the rail launch phase, RocketPy calculates reaction forces and internal bending moments at the rail button attachment points: + +**Rail Button Forces (N):** + +- ``rail_button1_normal_force`` : Normal reaction force at upper rail button +- ``rail_button1_shear_force`` : Shear (tangential) reaction force at upper rail button +- ``rail_button2_normal_force`` : Normal reaction force at lower rail button +- ``rail_button2_shear_force`` : Shear (tangential) reaction force at lower rail button + +**Rail Button Bending Moments (N⋅m):** + +- ``rail_button1_bending_moment`` : Time-dependent bending moment at upper rail button attachment +- ``max_rail_button1_bending_moment`` : Maximum absolute bending moment at upper rail button +- ``rail_button2_bending_moment`` : Time-dependent bending moment at lower rail button attachment +- ``max_rail_button2_bending_moment`` : Maximum absolute bending moment at lower rail button + +**Calculation Method:** + +Bending moments are calculated using beam theory assuming simple supports (rail buttons provide reaction forces but no moment reaction at rail contact). The total moment combines: + +1. Shear force × button height (cantilever moment from button standoff) +2. Normal force × distance to center of dry mass (lever arm effect) + +Moments are zero after rail departure and represent internal structural loads for airframe and fastener stress analysis. Requires ``button_height`` to be defined when adding rail buttons via ``rocket.set_rail_buttons()``. + +.. note:: + See Issue #893 for implementation details and validation approach. + + + Attitude and Orientation ~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index b2d3ae8a3..fc495ef3c 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -395,6 +395,70 @@ def angular_kinematics_data(self, *, filename=None): # pylint: disable=too-many plt.subplots_adjust(hspace=0.5) show_or_save_plot(filename) + def rail_buttons_bending_moments(self, *, filename=None): + """Prints out Rail Buttons Bending Moments graphs. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. Supported file endings are: + eps, jpg, jpeg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff + and webp (these are the formats supported by matplotlib). + + Returns + ------- + None + """ + if len(self.flight.rocket.rail_buttons) == 0: + print( + "No rail buttons were defined. Skipping rail button bending moment plots." + ) + elif self.flight.out_of_rail_time_index == 0: + print("No rail phase was found. Skipping rail button bending moment plots.") + else: + # Check if button_height is defined + rail_buttons_tuple = self.flight.rocket.rail_buttons[0] + if rail_buttons_tuple.component.button_height is None: + print("Rail button height not defined. Skipping bending moment plots.") + else: + plt.figure(figsize=(9, 3)) + + ax1 = plt.subplot(111) + ax1.plot( + self.flight.rail_button1_bending_moment[ + : self.flight.out_of_rail_time_index, 0 + ], + self.flight.rail_button1_bending_moment[ + : self.flight.out_of_rail_time_index, 1 + ], + label="Upper Rail Button", + ) + ax1.plot( + self.flight.rail_button2_bending_moment[ + : self.flight.out_of_rail_time_index, 0 + ], + self.flight.rail_button2_bending_moment[ + : self.flight.out_of_rail_time_index, 1 + ], + label="Lower Rail Button", + ) + ax1.set_xlim( + 0, + ( + self.flight.out_of_rail_time + if self.flight.out_of_rail_time > 0 + else self.flight.tFinal + ), + ) + ax1.legend() + ax1.grid(True) + ax1.set_xlabel("Time (s)") + ax1.set_ylabel("Bending Moment (N·m)") + ax1.set_title("Rail Button Bending Moments") + + show_or_save_plot(filename) + def rail_buttons_forces(self, *, filename=None): # pylint: disable=too-many-statements """Prints out all Rail Buttons Forces graphs available about the Flight. @@ -959,6 +1023,9 @@ def all(self): # pylint: disable=too-many-statements print("\n\nAerodynamic Forces Plots\n") self.aerodynamic_forces() + print("\n\nRail Buttons Bending Moments Plots\n") + self.rail_buttons_bending_moments() + print("\n\nRail Buttons Forces Plots\n") self.rail_buttons_forces() diff --git a/rocketpy/prints/flight_prints.py b/rocketpy/prints/flight_prints.py index d098fd776..1c0e3cb0c 100644 --- a/rocketpy/prints/flight_prints.py +++ b/rocketpy/prints/flight_prints.py @@ -358,6 +358,34 @@ def maximum_values(self): f"{self.flight.max_rail_button2_shear_force:.3f} N" ) + def rail_button_bending_moments(self): + """Prints rail button bending moment data. + + Returns + ------- + None + """ + if ( + len(self.flight.rocket.rail_buttons) == 0 + or self.flight.out_of_rail_time_index == 0 + ): + return + + # Check if button_height is defined + rail_buttons_tuple = self.flight.rocket.rail_buttons[0] + if rail_buttons_tuple.component.button_height is None: + return + + print("\nRail Button Bending Moments\n") + print( + "Maximum Upper Rail Button Bending Moment: " + f"{self.flight.max_rail_button1_bending_moment:.3f} N·m" + ) + print( + "Maximum Lower Rail Button Bending Moment: " + f"{self.flight.max_rail_button2_bending_moment:.3f} N·m" + ) + def stability_margin(self): """Prints out the stability margins of the flight at different times. @@ -429,5 +457,8 @@ def all(self): self.maximum_values() print() + self.rail_button_bending_moments() + print() + self.numerical_integration_settings() print() diff --git a/rocketpy/rocket/aero_surface/rail_buttons.py b/rocketpy/rocket/aero_surface/rail_buttons.py index c27e2949a..89331c99f 100644 --- a/rocketpy/rocket/aero_surface/rail_buttons.py +++ b/rocketpy/rocket/aero_surface/rail_buttons.py @@ -19,12 +19,19 @@ class RailButtons(AeroSurface): relative to one of the other principal axis. RailButtons.angular_position_rad : float Angular position of the rail buttons in radians. + RailButtons.button_height : float, optional + Height (standoff distance) of the rail button from the rocket + body surface to the rail contact point, in meters. Used for + calculating bending moments at the attachment point. + Default is None. If not provided, bending moments cannot be + calculated but flight dynamics remain unaffected. """ def __init__( self, buttons_distance, angular_position=45, + button_height=None, name="Rail Buttons", rocket_radius=None, ): @@ -48,6 +55,7 @@ def __init__( super().__init__(name, None, None) self.buttons_distance = buttons_distance self.angular_position = angular_position + self.button_height = button_height self.name = name self.rocket_radius = rocket_radius self.evaluate_lift_coefficient() @@ -104,6 +112,7 @@ def to_dict(self, **kwargs): # pylint: disable=unused-argument return { "buttons_distance": self.buttons_distance, "angular_position": self.angular_position, + "button_height": self.button_height, "name": self.name, "rocket_radius": self.rocket_radius, } @@ -113,6 +122,7 @@ def from_dict(cls, data): return cls( data["buttons_distance"], data["angular_position"], + data.get("button_height", None), data["name"], data["rocket_radius"], ) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 30ea66466..daa8719f8 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -470,6 +470,18 @@ class Flight: array. Flight.simulation_mode : str Simulation mode for the flight. Can be "6 DOF" or "3 DOF". + Flight.rail_button1_bending_moment : Function + Internal bending moment at upper rail button attachment point in N·m + as a function of time. Calculated using beam theory during rail phase. + Flight.max_rail_button1_bending_moment : float + Maximum internal bending moment experienced at upper rail button + attachment point during rail flight phase in N·m. + Flight.rail_button2_bending_moment : Function + Internal bending moment at lower rail button attachment point in N·m + as a function of time. Calculated using beam theory during rail phase. + Flight.max_rail_button2_bending_moment : float + Maximum internal bending moment experienced at lower rail button + attachment point during rail flight phase in N·m. """ def __init__( # pylint: disable=too-many-arguments,too-many-statements @@ -3958,3 +3970,142 @@ def __lt__(self, other): otherwise. """ return self.t < other.t + + @cached_property + def calculate_rail_button_bending_moments(self): + """ + Calculate internal bending moments at rail button attachment points. + + Uses beam theory to determine internal structural moments for stress + analysis of the rail button attachments (fasteners and airframe). + + The bending moment at each button attachment consists of: + 1. Bending from shear force at button contact point: M = S × h + where S is the shear (tangential) force and h is button height + 2. Direct moment contribution from the button's reaction forces + + Assumptions + ----------- + - Rail buttons act as simple supports: provide reaction forces (normal + and shear) but no moment reaction at the rail contact point. + - The rocket acts as a beam supported at two points (rail buttons). + - Bending moments arise from the lever arm effect of reaction forces + and the cantilever moment from button standoff height. + + The bending moment at each button attachment consists of: + 1. Normal force moment: M = N x d, where N is normal reaction force + and d is distance from button to center of dry mass + 2. Shear force cantilever moment: M = S x h, where S is shear force + and h is button standoff height + + Notes + ----- + - Calculated only during the rail phase of flight + - Maximum values use absolute values for worst-case stress analysis + - The bending moments represent internal stresses in the rocket + airframe at the rail button attachment points + + Returns + ------- + tuple + (rail_button1_bending_moment : Function, + max_rail_button1_bending_moment : float, + rail_button2_bending_moment : Function, + max_rail_button2_bending_moment : float) + + Where rail_button1/2_bending_moment are Function objects of time + in N·m, and max values are floats in N·m. + """ + # Check if rail buttons exist + null_moment = Function(0) + if len(self.rocket.rail_buttons) == 0: + warnings.warn( + "Trying to calculate rail button bending moments without " + "rail buttons defined. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + + # Get rail button geometry + rail_buttons_tuple = self.rocket.rail_buttons[0] + # Rail button standoff height + h_button = rail_buttons_tuple.component.button_height + if h_button is None: + warnings.warn( + "Rail button height not defined. Bending moments cannot be " + "calculated. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + upper_button_position = ( + rail_buttons_tuple.component.buttons_distance + + rail_buttons_tuple.position.z + ) + lower_button_position = rail_buttons_tuple.position.z + + # Get center of dry mass (handle both callable and property) + if callable(self.rocket.center_of_dry_mass_position): + cdm = self.rocket.center_of_dry_mass_position(self.rocket._csys) + else: + cdm = self.rocket.center_of_dry_mass_position + + # Distances from buttons to center of dry mass + d1 = abs(upper_button_position - cdm) + d2 = abs(lower_button_position - cdm) + + # forces + N1 = self.rail_button1_normal_force + N2 = self.rail_button2_normal_force + S1 = self.rail_button1_shear_force + S2 = self.rail_button2_shear_force + t = N1.source[:, 0] + + # Calculate bending moments at attachment points + # Primary contribution from shear force acting at button height + # Secondary contribution from normal force creating moment about attachment + m1_values = N2.source[:, 1] * d2 + S1.source[:, 1] * h_button + m2_values = N1.source[:, 1] * d1 + S2.source[:, 1] * h_button + + rail_button1_bending_moment = Function( + np.column_stack([t, m1_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + rail_button2_bending_moment = Function( + np.column_stack([t, m2_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + + # Maximum bending moments (absolute value for stress calculations) + max_rail_button1_bending_moment = float(np.max(np.abs(m1_values))) + max_rail_button2_bending_moment = float(np.max(np.abs(m2_values))) + + return ( + rail_button1_bending_moment, + max_rail_button1_bending_moment, + rail_button2_bending_moment, + max_rail_button2_bending_moment, + ) + + @property + def rail_button1_bending_moment(self): + """Upper rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[0] + + @property + def max_rail_button1_bending_moment(self): + """Maximum upper rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[1] + + @property + def rail_button2_bending_moment(self): + """Lower rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[2] + + @property + def max_rail_button2_bending_moment(self): + """Maximum lower rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[3] diff --git a/tests/unit/test_rail_buttons_bending_moments.py b/tests/unit/test_rail_buttons_bending_moments.py new file mode 100644 index 000000000..b61720d3d --- /dev/null +++ b/tests/unit/test_rail_buttons_bending_moments.py @@ -0,0 +1,222 @@ +"""Unit tests for RailButtons bending moment formulas and calculation logic. + +These tests focus on verifying the mathematical correctness and realism +of the calculate_rail_button_bending_moments method in isolation. +""" + +import warnings + +import matplotlib.pyplot as plt +import numpy as np + +from rocketpy import Environment, Flight +from rocketpy.mathutils.function import Function +from rocketpy.rocket.aero_surface.rail_buttons import RailButtons + + +def test_bending_moment_zero_without_rail_buttons(calisto_motorless): + """Verify bending moments return zero functions if no rail buttons present. + + Parameters + ---------- + calisto_motorless : rocketpy.Rocket + Basic rocket without rail buttons. + """ + # Create a flight with this rocket (no rail buttons) + env = Environment(latitude=0, longitude=0) + env.set_atmospheric_model(type="standard_atmosphere") + flight = Flight(rocket=calisto_motorless, environment=env, rail_length=1) + + # Should return zero moments + moments = flight.calculate_rail_button_bending_moments + m1_func, max_m1, m2_func, max_m2 = moments + + # Verify types + assert isinstance(m1_func, Function) + assert isinstance(m2_func, Function) + assert isinstance(max_m1, float) + assert isinstance(max_m2, float) + + # Verify zero functions + assert m1_func(0) == 0 + assert m1_func(1) == 0 + assert m2_func(0) == 0 + assert m2_func(1) == 0 + assert max_m1 == 0.0 + assert max_m2 == 0.0 + + +def test_bending_moment_zero_with_none_button_height(calisto_motorless): + """Verify bending moments return zero when button_height is explicitly None. + + Parameters + ---------- + calisto_motorless : rocketpy.Rocket + Basic rocket fixture. + """ + + # Create rail buttons, then explicitly set height to None + calisto_motorless.set_rail_buttons( + upper_button_position=0.5, + lower_button_position=-0.5, + angular_position=45, + ) + calisto_motorless.rail_buttons[0].component.button_height = None # explicit None + + # Sanity check + assert calisto_motorless.rail_buttons[0].component.button_height is None + + env = Environment(latitude=0, longitude=0) + env.set_atmospheric_model(type="standard_atmosphere") + flight = Flight(rocket=calisto_motorless, environment=env, rail_length=1) + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + moments = flight.calculate_rail_button_bending_moments + assert len(w) == 1 + assert "button height not defined" in str(w[0].message).lower() + + m1_func, max_m1, _, max_m2 = moments + assert m1_func(0) == 0 + assert max_m1 == 0.0 + assert max_m2 == 0.0 + + +def test_bending_moment_zero_with_default_button_height(calisto_motorless): + """Verify bending moments return zero when button_height uses default value. + + Parameters + ---------- + calisto_motorless : rocketpy.Rocket + Basic rocket fixture. + """ + # Add rail buttons without specifying button_height (tests default) + calisto_motorless.set_rail_buttons( + upper_button_position=0.5, + lower_button_position=-0.5, + angular_position=45, + # button_height not specified - should default to None + ) + + # Verify default is None + assert calisto_motorless.rail_buttons[0].component.button_height is None + + +def test_railbuttons_serialization_roundtrip(): + """Test RailButtons to_dict and from_dict serialization roundtrip.""" + rb_orig = RailButtons( + buttons_distance=0.7, + angular_position=60, + button_height=0.02, + name="Test Rail Buttons", + rocket_radius=0.045, + ) + + data = rb_orig.to_dict() + rb_reconstructed = RailButtons.from_dict(data) + + assert rb_reconstructed.buttons_distance == rb_orig.buttons_distance + assert rb_reconstructed.angular_position == rb_orig.angular_position + assert rb_reconstructed.button_height == rb_orig.button_height + assert rb_reconstructed.name == rb_orig.name + assert rb_reconstructed.rocket_radius == rb_orig.rocket_radius + + +def test_railbuttons_serialization_backward_compat(): + """Test backward compatibility when button_height is missing from dict.""" + # Simulate old serialized data without button_height + old_data = { + "buttons_distance": 0.5, + "angular_position": 45, + "name": "Legacy Buttons", + "rocket_radius": 0.05, + } + + rb = RailButtons.from_dict(old_data) + assert rb.button_height is None # Should use None default + + +def test_railbuttons_angular_position_rad_property(): + """Test angular_position_rad property calculation.""" + rb = RailButtons(buttons_distance=0.5, angular_position=30) + expected_rad = np.radians(30) + assert np.isclose(rb.angular_position_rad, expected_rad, rtol=1e-10) + + +def test_railbuttons_no_aero_contribution(): + """Test RailButtons provide zero aerodynamic contributions.""" + rb = RailButtons(buttons_distance=0.5) + + rb.evaluate_center_of_pressure() + assert rb.cp == (0, 0, 0) + + rb.evaluate_lift_coefficient() + assert rb.clalpha(1.0) == 0 # Zero lift derivative + assert rb.cl(0.1, 1.0) == 0 # Zero lift coefficient + + +def test_rail_button_bending_moments_prints(flight_calisto_robust, capsys): + """Test that bending moments are printed correctly in flight.prints. + + Parameters + ---------- + flight_calisto_robust : rocketpy.Flight + Flight fixture with motor and rail buttons. + capsys : pytest fixture + Captures stdout/stderr. + """ + # Set button height on existing rail buttons + flight_calisto_robust.rocket.rail_buttons[0].component.button_height = 0.02 + + # Call the print method + flight_calisto_robust.prints.rail_button_bending_moments() + + # Capture output + captured = capsys.readouterr() + + # Verify output contains bending moment data + assert "Rail Button Bending Moments" in captured.out + assert "Maximum Upper Rail Button Bending Moment" in captured.out + assert "Maximum Lower Rail Button Bending Moment" in captured.out + assert "N·m" in captured.out + + +def test_rail_button_bending_moments_plot_with_height(flight_calisto_robust): + """Test that bending moments plot is created when button_height is defined. + + Parameters + ---------- + flight_calisto_robust : rocketpy.Flight + Flight fixture with motor and rail buttons. + """ + # Set button height on existing rail buttons + flight_calisto_robust.rocket.rail_buttons[0].component.button_height = 0.02 + + # Should not raise an error + flight_calisto_robust.plots.rail_buttons_bending_moments() + + # Close the figure to avoid memory issues + plt.close("all") + + +def test_rail_button_bending_moments_plot_without_height(flight_calisto_robust, capsys): + """Test that bending moments plot is skipped when button_height is None. + + Parameters + ---------- + flight_calisto_robust : rocketpy.Flight + Flight fixture with motor and rail buttons. + capsys : pytest fixture + Captures stdout/stderr. + """ + # Ensure button_height is None (it should be by default now) + flight_calisto_robust.rocket.rail_buttons[0].component.button_height = None + + # Call plot method + flight_calisto_robust.plots.rail_buttons_bending_moments() + + # Capture output + captured = capsys.readouterr() + + # Should print skip message + assert "height not defined" in captured.out