ENH: Parachute Inflation Modeling #867
Conversation
|
@ArthurJWH please rebase your branch to develop so we can properly review this PR. I'm under the impression the 26 commits of this PR should not belong to this git history |
There was a problem hiding this comment.
Pull Request Overview
This PR improves the parachute model by adding geometric parametrization to represent parachutes as hemispheroids with radius, height, and porosity properties. The key changes include:
- Added
radius,height, andporosityparameters to theParachuteclass with corresponding added-mass coefficient calculation - Updated the flight simulation's parachute dynamics to include inflation modeling and added-mass calculation based on geometric properties
- Extended
StochasticParachuteto support randomization of the new geometric parameters - Updated test assertions to reflect the changes in physics calculations
Reviewed Changes
Copilot reviewed 10 out of 10 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| rocketpy/rocket/parachute.py | Adds radius, height, porosity attributes and calculates added-mass coefficient from porosity |
| rocketpy/rocket/rocket.py | Updates add_parachute method to accept new geometric parameters |
| rocketpy/simulation/flight.py | Implements inflation modeling and added-mass calculation using parachute geometry; updates parachute callbacks |
| rocketpy/stochastic/stochastic_parachute.py | Extends stochastic model to support randomization of radius, height, and porosity |
| tests/unit/test_utilities.py | Updates safety factor assertion to match new physics calculations |
| tests/unit/test_flight.py | Updates aerodynamic forces test values to reflect changed dynamics |
| docs/user/rocket/rocket_usage.rst | Adds example usage of new parachute parameters |
| docs/user/first_simulation.rst | Adds example usage of new parachute parameters |
| README.md | Adds example usage of new parachute parameters |
| CHANGELOG.md | Documents the enhancement |
ArthurJWH
force-pushed
the
enh/parachute-inflation-modeling
branch
from
e02aceb to
8cef755
Compare
Gui-FernandesBR
force-pushed
the
enh/parachute-inflation-modeling
branch
from
3da56ff to
b861b2d
Compare
|
@ArthurJWH tests and linters are not passing.... and, finally, I think there's smth wrong with your git tree again... Could you check the git diff and verify everything is good with your PR? Let us know whenever you have the PR ready to be reviewed again please. |
…/RocketPy-Team/RocketPy into enh/parachute-inflation-modeling
| ), | ||
| lambda self: setattr( | ||
| self, | ||
| "parachute_volume", |
There was a problem hiding this comment.
@ArthurJWH why do you need to calculate the parachute_volume within the flight simulation? Can you calculate and store it in the Parachute class?
There was a problem hiding this comment.
My original idea was to assume the initial radius being the rocket radius (compartment from where it is ejected), but I cannot access the rocket radius within Parachute class. I can add a new parameter "initial radius" for the Parachute class, if you think that is better
| free_stream_speed = (freestream_x**2 + freestream_y**2 + freestream_z**2) ** 0.5 | ||
|
|
||
| # Initialize parachute geometrical parameters | ||
| radius = self.parachute_radius |
There was a problem hiding this comment.
Wha is the benefit behind this radius = self.parachute_radius?
There was a problem hiding this comment.
I thought it would help with the readability of the code, especially when the parameters are used in the convoluted equations. I can remove it
There was a problem hiding this comment.
I understand it.
it does not contribute so much to the readability actually.
We can save a few 0.0001 ns by removing this line
| # Initialize parachute geometrical parameters | ||
| radius = self.parachute_radius | ||
| height = self.parachute_height | ||
| inflated_radius = ( |
There was a problem hiding this comment.
calculating inflated_radius in the Flight class? Do you really need it?
There was a problem hiding this comment.
inflated_radius is the parachute radius in the current Flight stage, it is used to calculate the surface area of the parachute
| ) | ||
| ) | ||
|
|
||
| # Getting time step |
There was a problem hiding this comment.
do you really need to store step.t0? Why?
Should this really be a public attribute?
There was a problem hiding this comment.
I used it to have access to the time step between u_dot_parachute calls (not ideal). I can change it for a private attribute, or do you have any other suggestion?
| volume_flow = ( | ||
| rho | ||
| * freestream_z # considering parachute as vertical | ||
| * ( | ||
| (math.pi * inflated_radius**2) | ||
| - (self.parachute_porosity * surface_area) | ||
| ) | ||
| ) |
There was a problem hiding this comment.
The volume_flow calculation has incorrect units and logic. Currently, it computes rho * velocity * area, which yields mass flow rate (kg/s), not volume flow rate (m³/s). The variable should be renamed to 'mass_flow' and the integration on line 2779 should add mass, not volume. Alternatively, remove the rho multiplication here to get true volume flow, then multiply by rho only when computing the added mass on line 2782.
| e = math.sqrt(1 - (inflated_height**2) / (inflated_radius**2)) | ||
| surface_area = ( | ||
| math.pi * inflated_radius**2 * (1 + (1 - e**2) / e * math.atanh(e)) | ||
| ) | ||
| else: | ||
| e = math.sqrt(1 - (inflated_radius**2) / (inflated_height**2)) | ||
| surface_area = ( | ||
| math.pi | ||
| * inflated_radius**2 | ||
| * (1 + inflated_height / (e * inflated_radius) * math.asin(e)) |
There was a problem hiding this comment.
The surface area calculation can fail with a domain error when the eccentricity formula yields a negative value under sqrt. This can occur if inflated_height > inflated_radius when radius > height (line 2751), or if inflated_radius > inflated_height when radius <= height (line 2756). Add validation to ensure the discriminant is non-negative before taking the square root, or clamp the eccentricity to valid ranges [0, 1).
| e = math.sqrt(1 - (inflated_height**2) / (inflated_radius**2)) | |
| surface_area = ( | |
| math.pi * inflated_radius**2 * (1 + (1 - e**2) / e * math.atanh(e)) | |
| ) | |
| else: | |
| e = math.sqrt(1 - (inflated_radius**2) / (inflated_height**2)) | |
| surface_area = ( | |
| math.pi | |
| * inflated_radius**2 | |
| * (1 + inflated_height / (e * inflated_radius) * math.asin(e)) | |
| discriminant = 1 - (inflated_height**2) / (inflated_radius**2) | |
| if discriminant < 0: | |
| warnings.warn( | |
| f"Parachute eccentricity discriminant negative ({discriminant:.6f}); " | |
| "clamping to zero. Check parachute geometry parameters." | |
| ) | |
| discriminant = 0.0 | |
| elif discriminant > 1: | |
| warnings.warn( | |
| f"Parachute eccentricity discriminant > 1 ({discriminant:.6f}); " | |
| "clamping to one. Check parachute geometry parameters." | |
| ) | |
| discriminant = 1.0 | |
| e = math.sqrt(discriminant) | |
| surface_area = ( | |
| math.pi * inflated_radius**2 * (1 + (1 - e**2) / e * math.atanh(e)) if e != 0 else | |
| math.pi * inflated_radius**2 * 2 # limit as e->0 | |
| ) | |
| else: | |
| discriminant = 1 - (inflated_radius**2) / (inflated_height**2) | |
| if discriminant < 0: | |
| warnings.warn( | |
| f"Parachute eccentricity discriminant negative ({discriminant:.6f}); " | |
| "clamping to zero. Check parachute geometry parameters." | |
| ) | |
| discriminant = 0.0 | |
| elif discriminant > 1: | |
| warnings.warn( | |
| f"Parachute eccentricity discriminant > 1 ({discriminant:.6f}); " | |
| "clamping to one. Check parachute geometry parameters." | |
| ) | |
| discriminant = 1.0 | |
| e = math.sqrt(discriminant) | |
| surface_area = ( | |
| math.pi | |
| * inflated_radius**2 | |
| * (1 + inflated_height / (e * inflated_radius) * math.asin(e)) if e != 0 else | |
| math.pi * inflated_radius**2 * 2 # limit as e->0 |
There was a problem hiding this comment.
Since the inflated radius / inflated height = radius / height, the if else statement should already account for sign of the value under the square root
| # Getting time step | ||
| self.t0 = getattr(self, "t0", t) | ||
| t1 = t | ||
| dt = t1 - self.t0 | ||
|
|
||
| # Integrating parachute volume | ||
| self.parachute_volume += volume_flow * dt | ||
|
|
||
| # Dragged air mass | ||
| ma = self.parachute_volume * rho | ||
|
|
||
| # Moving time step | ||
| self.t0 = t1 |
There was a problem hiding this comment.
This manual integration approach using timestep differences is inefficient and can introduce numerical errors. The ODE solver will call this derivative function multiple times per timestep with different intermediate values of t, but the integration assumes each call represents a true timestep. This will cause parachute_volume to accumulate incorrectly. Consider reformulating this as a proper differential equation by adding parachute_volume to the state vector u, and returning its time derivative as part of u_dot.
| # Calculate volume flow of air into parachute | ||
| volume_flow = ( | ||
| rho | ||
| * freestream_z # considering parachute as vertical |
There was a problem hiding this comment.
The volume flow calculation assumes the parachute is vertical by using only freestream_z. However, parachutes don't necessarily descend perfectly vertically, especially in windy conditions or during initial deployment. The horizontal velocity components (freestream_x, freestream_y) should also contribute to the flow into the parachute. Consider using the total freestream velocity magnitude or projecting the velocity onto the parachute's orientation axis.
| # Calculate volume flow of air into parachute | |
| volume_flow = ( | |
| rho | |
| * freestream_z # considering parachute as vertical | |
| # Calculate volume flow of air into parachute using total freestream velocity magnitude | |
| freestream_velocity_magnitude = math.sqrt( | |
| freestream_x**2 + freestream_y**2 + freestream_z**2 | |
| ) | |
| volume_flow = ( | |
| rho | |
| * freestream_velocity_magnitude |
| (4 / 3) | ||
| * math.pi | ||
| * (self.parachute_height / self.parachute_radius) | ||
| * ( | ||
| min( | ||
| self.parachute_radius, | ||
| self.rocket.radius, | ||
| ) | ||
| ) | ||
| * 3, |
There was a problem hiding this comment.
The parachute_volume initialization formula appears to be incorrect. The expression simplifies to 4 * π * (height/radius) * min(radius, rocket.radius), which has units of m² (area), not m³ (volume). For a hemispheroid parachute, the initial volume should be approximately (2/3) * π * radius^2 * height (half of an oblate spheroid). The current formula will produce incorrect initial volume values that will affect the entire inflation simulation.
| (4 / 3) | |
| * math.pi | |
| * (self.parachute_height / self.parachute_radius) | |
| * ( | |
| min( | |
| self.parachute_radius, | |
| self.rocket.radius, | |
| ) | |
| ) | |
| * 3, | |
| (2 / 3) | |
| * math.pi | |
| * self.parachute_radius ** 2 | |
| * self.parachute_height, |
| ) | ||
| * 3, | ||
| ), | ||
| lambda self: delattr(self, "t0"), |
There was a problem hiding this comment.
Attempting to delete the 't0' attribute before it has been set will raise an AttributeError. The 't0' attribute is only created later in u_dot_parachute (line 2774) using getattr with a default value. This callback should either be removed or moved to execute after the parachute phase ends, not when it begins.
| lambda self: delattr(self, "t0"), | |
| lambda self: delattr(self, "t0") if hasattr(self, "t0") else None, |
| (4 / 3) | ||
| * math.pi | ||
| * (self.parachute_height / self.parachute_radius) | ||
| * ( | ||
| min( | ||
| self.parachute_radius, | ||
| self.rocket.radius, | ||
| ) | ||
| ) | ||
| * 3, |
There was a problem hiding this comment.
The parachute_volume initialization formula appears to be incorrect. The expression simplifies to 4 * π * (height/radius) * min(radius, rocket.radius), which has units of m² (area), not m³ (volume). For a hemispheroid parachute, the initial volume should be approximately (2/3) * π * radius^2 * height (half of an oblate spheroid). The current formula will produce incorrect initial volume values that will affect the entire inflation simulation.
| (4 / 3) | |
| * math.pi | |
| * (self.parachute_height / self.parachute_radius) | |
| * ( | |
| min( | |
| self.parachute_radius, | |
| self.rocket.radius, | |
| ) | |
| ) | |
| * 3, | |
| (2 / 3) | |
| * math.pi | |
| * (min(self.parachute_radius, self.rocket.radius) ** 2) | |
| * self.parachute_height, |
| lambda self: setattr( | ||
| self, | ||
| "parachute_volume", | ||
| (4 / 3) | ||
| * math.pi | ||
| * (self.parachute_height / self.parachute_radius) | ||
| * ( | ||
| min( | ||
| self.parachute_radius, | ||
| self.rocket.radius, | ||
| ) | ||
| ) | ||
| * 3, | ||
| ), | ||
| lambda self: delattr(self, "t0"), |
There was a problem hiding this comment.
This code block is duplicated at lines 768-782. Consider extracting the parachute initialization callbacks into a helper method to avoid duplication and ensure both locations remain consistent. The duplication increases maintenance burden and the risk of inconsistencies if one location is updated but not the other.
| # Initialize parachute geometrical parameters | ||
| radius = self.parachute_radius | ||
| height = self.parachute_height | ||
| inflated_radius = ( | ||
| (3 * self.parachute_volume * radius) / (4 * math.pi * height) | ||
| ) ** (1 / 3) | ||
| inflated_height = inflated_radius * height / radius | ||
|
|
||
| # Calculate the surface area of the parachute | ||
| if radius > height: | ||
| e = math.sqrt(1 - (inflated_height**2) / (inflated_radius**2)) | ||
| surface_area = ( | ||
| math.pi * inflated_radius**2 * (1 + (1 - e**2) / e * math.atanh(e)) | ||
| ) | ||
| else: | ||
| e = math.sqrt(1 - (inflated_radius**2) / (inflated_height**2)) | ||
| surface_area = ( | ||
| math.pi | ||
| * inflated_radius**2 | ||
| * (1 + inflated_height / (e * inflated_radius) * math.asin(e)) | ||
| ) | ||
|
|
||
| # Calculate volume flow of air into parachute | ||
| volume_flow = ( | ||
| rho | ||
| * freestream_z # considering parachute as vertical | ||
| * ( | ||
| (math.pi * inflated_radius**2) | ||
| - (self.parachute_porosity * surface_area) | ||
| ) | ||
| ) | ||
|
|
||
| # Getting time step | ||
| self.t0 = getattr(self, "t0", t) | ||
| t1 = t | ||
| dt = t1 - self.t0 | ||
|
|
||
| # Integrating parachute volume | ||
| self.parachute_volume += volume_flow * dt | ||
|
|
||
| # Dragged air mass | ||
| ma = self.parachute_volume * rho | ||
|
|
||
| # Moving time step | ||
| self.t0 = t1 |
There was a problem hiding this comment.
The new parachute inflation model (lines 2741-2785) lacks test coverage. This is a significant new feature that modifies how parachute dynamics are simulated. Consider adding tests that verify: (1) parachute_volume increases over time during inflation, (2) the volume converges to the expected inflated volume, (3) the surface area calculations produce physically reasonable values, and (4) edge cases like very small or very large parachutes are handled correctly.
| ) | ||
| * 3, | ||
| ), | ||
| lambda self: delattr(self, "t0"), |
There was a problem hiding this comment.
Attempting to delete the 't0' attribute before it has been set will raise an AttributeError. The 't0' attribute is only created later in u_dot_parachute (line 2774) using getattr with a default value. This callback should either be removed or moved to execute after the parachute phase ends, not when it begins.
| lambda self: delattr(self, "t0"), | |
| lambda self: delattr(self, "t0") if hasattr(self, "t0") else None, |
Pull request type
Checklist
black rocketpy/ tests/) has passed locallypytest tests -m slow --runslow) have passed locallyCHANGELOG.mdhas been updated (if relevant)Current behavior
When triggered, the parachute is modeled as fully opened
New behavior
The parachute inflation will be simulated using a volumetric flow model with simple inviscid flow
Breaking change