Skip to content

Commit

Permalink
MNT: improve overall evaluate_geom
Browse files Browse the repository at this point in the history
  • Loading branch information
@MateusStano
MateusStano committed Sep 17, 2024
1 parent 3845316 commit 2869e91
Showing 1 changed file with 84 additions and 78 deletions.
162 changes: 84 additions & 78 deletions rocketpy/rocket/aero_surface/fins/free_form_fins.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,87 +197,91 @@ def evaluate_center_of_pressure(self):
self.cp = (self.cpx, self.cpy, self.cpz)

def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements
"""Calculates and saves fin set's geometrical parameters such as the
fins' area, aspect ratio and parameters for roll movement. The
calculations related to the free form fins points are done in the exact
same way as Open Rocket.
"""
Calculates and saves the fin set's geometrical parameters such as the
fin area, aspect ratio, and parameters related to roll movement. This method
uses similar calculations to those in OpenRocket for free-form fin shapes.
Returns
-------
None
"""
# pylint: disable=invalid-name
# Calculate the fin area (Af) using the Shoelace theorem (polygon area formula)
Af = 0
for i in range(len(self.shape_points) - 1):
Af += (self.shape_points[i][1] + self.shape_points[i + 1][1]) * (
self.shape_points[i][0] - self.shape_points[i + 1][0]
)
x1, y1 = self.shape_points[i]
x2, y2 = self.shape_points[i + 1]
Af += (y1 + y2) * (x1 - x2)
Af = abs(Af) / 2
if Af < 1e-6:
raise ValueError("Fin area is too small. Check the shape_points.")

AR = 2 * self.span**2 / Af # Fin aspect ratio
# Calculate aspect ratio (AR) and lift interference factors
AR = 2 * self.span**2 / Af # Aspect ratio
tau = (self.span + self.rocket_radius) / self.rocket_radius
lift_interference_factor = 1 + 1 / tau

# Calculate roll forcing interference factor using OpenRocket's approach
roll_forcing_interference_factor = (1 / np.pi**2) * (
(np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2)
+ ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2))
* np.arcsin((tau**2 - 1) / (tau**2 + 1))
- (2 * np.pi * (tau + 1)) / (tau * (tau - 1))
+ ((tau**2 + 1) ** 2)
/ (tau**2 * (tau - 1) ** 2)
+ ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2))
* (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2
- (4 * (tau + 1))
/ (tau * (tau - 1))
* np.arcsin((tau**2 - 1) / (tau**2 + 1))
+ (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau))
)

# Discretize fin shape interpolation into points_per_line points
points_per_line = 40 # Arbitrary number of points per line, same as OPR
# Define number of interpolation points along the span of the fin
points_per_line = 40 # Same as OpenRocket

chord_lead = np.ones(points_per_line) * np.inf # Leading edge of the chord
chord_trail = np.ones(points_per_line) * -np.inf # Trailing edge of the chord
chord_length = np.zeros(points_per_line) # Chord length
# Initialize arrays for leading/trailing edge and chord lengths
chord_lead = np.ones(points_per_line) * np.inf # Leading edge x coordinates
chord_trail = np.ones(points_per_line) * -np.inf # Trailing edge x coordinates
chord_length = np.zeros(
points_per_line
) # Chord length for each spanwise section

# Calculate chord length and leading/trailing edge
# If fin is jagged, calculation will be for the equivalent shape
# Discretize fin shape and calculate chord length, leading, and trailing edges
for p in range(1, len(self.shape_points)):
# Coordinates of the two points
x1 = np.float64(self.shape_points[p - 1][0])
y1 = np.float64(self.shape_points[p - 1][1])
x2 = np.float64(self.shape_points[p][0])
y2 = np.float64(self.shape_points[p][1])

# Correction for jagged fin
previous_point = int(y1 * 1.0001 / self.span * (points_per_line - 1))
current_point = int(y2 * 1.0001 / self.span * (points_per_line - 1))
previous_point = max(min(previous_point, points_per_line - 1), 0)
current_point = max(min(current_point, points_per_line - 1), 0)
if previous_point > current_point:
previous_point, current_point = current_point, previous_point

# Calculate chord length and leading/trailing edge
for i in range(previous_point, current_point + 1):
x1, y1 = self.shape_points[p - 1]
x2, y2 = self.shape_points[p]

# Compute corresponding points along the fin span (clamp to valid range)
prev_idx = int(y1 * 1.0001 / self.span * (points_per_line - 1))
curr_idx = int(y2 * 1.0001 / self.span * (points_per_line - 1))
prev_idx = np.clip(prev_idx, 0, points_per_line - 1)
curr_idx = np.clip(curr_idx, 0, points_per_line - 1)

if prev_idx > curr_idx:
prev_idx, curr_idx = curr_idx, prev_idx

# Compute intersection of fin edge with each spanwise section
for i in range(prev_idx, curr_idx + 1):
y = i * self.span / (points_per_line - 1)
x = max(
min(
if y1 != y2:
x = np.clip(
(y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2,
min(x1, x2),
max(x1, x2),
),
min(x1, x2),
)
if x < chord_lead[i]:
chord_lead[i] = x
if x > chord_trail[i]:
chord_trail[i] = x
)
else:
x = x1 # Handle horizontal segments

# Update leading and trailing edge positions
chord_lead[i] = min(chord_lead[i], x)
chord_trail[i] = max(chord_trail[i], x)

# Update chord length
if y1 < y2:
chord_length[i] -= x
else:
chord_length[i] += x

# Remove infinities and nans for cases where there are equal coordinates
# Replace infinities and handle invalid values in chord_lead and chord_trail
for i in range(points_per_line):
if (
np.isinf(chord_lead[i])
Expand All @@ -292,67 +296,69 @@ def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statement
if chord_length[i] > chord_trail[i] - chord_lead[i]:
chord_length[i] = chord_trail[i] - chord_lead[i]

# Initialize integration variables
radius = self.rocket_radius # Rocket radius
area = 0 # Area of the fin, this will be different from Af if fin is jagged
mac_length = 0 # Mean aerodynamic chord length)
mac_lead = 0 # Mean aerodynamic chord leading edge x coordinate
mac_span = 0 # Mean aerodynamic chord span wise position (Yma)
cos_gamma = 0 # Cosine of the sweep angle
# Initialize integration variables for various aerodynamic and roll properties
radius = self.rocket_radius
total_area = 0
mac_length = 0 # Mean aerodynamic chord length
mac_lead = 0 # Mean aerodynamic chord leading edge
mac_span = 0 # Mean aerodynamic chord spanwise position (Yma)
cos_gamma_sum = 0 # Sum of cosine of the sweep angle
roll_geometrical_constant = 0
roll_damping_interference_factor_numerator = 0
roll_damping_interference_factor_denominator = 0
roll_damping_numerator = 0
roll_damping_denominator = 0

# Perform integration
# Perform integration over spanwise sections
dy = self.span / (points_per_line - 1)
for i in range(points_per_line):
length = chord_trail[i] - chord_lead[i]
chord = chord_trail[i] - chord_lead[i]
y = i * dy

mac_length += length * length
mac_span += y * length
mac_lead += chord_lead[i] * length
area += length
# Update integration variables
mac_length += chord * chord
mac_span += y * chord
mac_lead += chord_lead[i] * chord
total_area += chord
roll_geometrical_constant += chord_length[i] * (radius + y) ** 2
roll_damping_interference_factor_numerator += (
radius**3 * length / (radius + y) ** 2
)
roll_damping_interference_factor_denominator += (radius + y) * length
roll_damping_numerator += radius**3 * chord / (radius + y) ** 2
roll_damping_denominator += (radius + y) * chord

# Update cosine of sweep angle (cos_gamma)
if i > 0:
dx = (chord_trail[i] + chord_lead[i]) / 2 - (
chord_trail[i - 1] + chord_lead[i - 1]
) / 2
cos_gamma += dy / (dx**2 + dy**2) ** 0.5
cos_gamma_sum += dy / np.hypot(dx, dy)

# Finalize mean aerodynamic chord properties
mac_length *= dy
mac_span *= dy
mac_lead *= dy
area *= dy
total_area *= dy
roll_geometrical_constant *= dy
roll_damping_interference_factor_numerator *= dy
roll_damping_interference_factor_denominator *= dy
mac_length /= area
mac_span /= area
mac_lead /= area
cos_gamma /= points_per_line - 1

# Store values
roll_damping_numerator *= dy
roll_damping_denominator *= dy

mac_length /= total_area
mac_span /= total_area
mac_lead /= total_area
cos_gamma = cos_gamma_sum / (points_per_line - 1)

# Store computed values
self.Af = Af # Fin area
self.AR = AR # Aspect Ratio
self.gamma_c = np.arccos(cos_gamma) # Mid chord angle
self.Yma = mac_span # Span wise coord of mean aero chord
self.AR = AR # Aspect ratio
self.gamma_c = np.arccos(cos_gamma) # Sweep angle
self.Yma = mac_span # Mean aerodynamic chord spanwise position
self.mac_length = mac_length
self.mac_lead = mac_lead
self.tau = tau
self.roll_geometrical_constant = roll_geometrical_constant
self.lift_interference_factor = lift_interference_factor
self.roll_forcing_interference_factor = roll_forcing_interference_factor
self.roll_damping_interference_factor = 1 + (
roll_damping_interference_factor_numerator
/ roll_damping_interference_factor_denominator
roll_damping_numerator / roll_damping_denominator
)

# Evaluate the shape and finalize geometry
self.evaluate_shape()

def evaluate_shape(self):
Expand Down

0 comments on commit 2869e91

Please sign in to comment.