From 2869e91399e0d170298818252c57f46bf13a48c0 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Wed, 18 Sep 2024 01:22:57 +0200 Subject: [PATCH] MNT: improve overall evaluate_geom --- .../aero_surface/fins/free_form_fins.py | 162 +++++++++--------- 1 file changed, 84 insertions(+), 78 deletions(-) diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py index 3aeb13d7e..0b24bf76e 100644 --- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py +++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py @@ -197,35 +197,37 @@ 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)) @@ -233,51 +235,53 @@ def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statement + (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]) @@ -292,56 +296,58 @@ 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 @@ -349,10 +355,10 @@ def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statement 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):