From acbfc2c3c06597c251777ae254db6a3a491bb593 Mon Sep 17 00:00:00 2001 From: jellepoland Date: Fri, 26 Sep 2025 09:50:27 +0200 Subject: [PATCH 1/3] changed tutorial --- ...etry.yaml => struc_geometry_surfplan.yaml} | 66 +++++++++---------- examples/TUDELFT_V3_KITE/tutorial.py | 9 ++- 2 files changed, 37 insertions(+), 38 deletions(-) rename data/TUDELFT_V3_KITE/CAD_derived_geometry/{struc_geometry.yaml => struc_geometry_surfplan.yaml} (94%) diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry_surfplan.yaml similarity index 94% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry_surfplan.yaml index 61acaac..86ea16e 100644 --- a/data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry.yaml +++ b/data/TUDELFT_V3_KITE/CAD_derived_geometry/struc_geometry_surfplan.yaml @@ -217,72 +217,72 @@ bridle_particles: - [87, 1.167853, 1.92072, 3.288303] - [88, 1.204082, 0.654966, 3.540403] bridle_connections: - headers: [name, ci, cj, ck] + headers: [name, ci, cj] data: - [a1, 21, 23] - - [a1, 35, 37] + - [a1, 55, 57] - [a2, 22, 25] - - [a2, 36, 39] + - [a2, 56, 59] - [a3, 26, 28] - - [a3, 40, 42] + - [a3, 60, 62] - [a4, 33, 37] - - [a4, 47, 51] + - [a4, 67, 71] - [A1, 23, 30] - - [A1, 37, 44] + - [A1, 57, 64] - [A2, 25, 31] - - [A2, 39, 45] + - [A2, 59, 65] - [A3, 28, 38] - - [A3, 42, 52] + - [A3, 62, 72] - [A4, 37, 36] - - [A4, 51, 50] + - [A4, 71, 70] - [AI, 30, 32] - - [AI, 44, 46] + - [AI, 64, 66] - [AII, 31, 32] - - [AII, 45, 46] + - [AII, 65, 66] - [AIII, 38, 36] - - [AIII, 52, 50] + - [AIII, 72, 70] - [amain, 32, 34] - - [amain, 46, 48] + - [amain, 66, 68] - [amain2, 36, 34] - - [amain2, 50, 48] + - [amain2, 70, 68] - [a5.1, 34, 35] - - [a5.1, 48, 49] + - [a5.1, 68, 69] - [b1, 24, 23] - - [b1, 38, 37] + - [b1, 58, 57] - [b2, 27, 25] - - [b2, 41, 39] + - [b2, 61, 59] - [b3, 29, 28] - - [b3, 43, 42] + - [b3, 63, 62] - [b4, 39, 37] - - [b4, 53, 51] + - [b4, 73, 71] - [e1, 40, 30] - - [e1, 54, 44] + - [e1, 74, 64] - [e2, 41, 31] - - [e2, 55, 45] + - [e2, 75, 65] - [e3, 42, 38] - - [e3, 56, 52] + - [e3, 76, 72] - [br1, 54, 47] - - [br1, 68, 61] + - [br1, 88, 81] - [br2, 53, 47] - - [br2, 67, 61] + - [br2, 87, 81] - [br3, 52, 46] - - [br3, 66, 60] + - [br3, 86, 80] - [br4, 51, 46] - - [br4, 65, 60] + - [br4, 85, 80] - [br5, 48, 49] - - [br5, 62, 63] + - [br5, 82, 83] - [br6, 50, 49] - - [br6, 64, 63] + - [br6, 84, 83] - [BR1, 47, 45] - - [BR1, 61, 59] + - [BR1, 81, 79] - [BR2, 46, 45] - - [BR2, 60, 59] + - [BR2, 80, 79] - [BR3, 49, 44] - - [BR3, 63, 58] + - [BR3, 83, 78] - [BRI, 45, 44] - - [BRI, 59, 58] + - [BRI, 79, 78] - [brmain, 44, 43] - - [brmain, 58, 57] + - [brmain, 78, 77] bridle_lines: headers: [name, rest_length, diameter, material, density] data: diff --git a/examples/TUDELFT_V3_KITE/tutorial.py b/examples/TUDELFT_V3_KITE/tutorial.py index 38a01b6..666241b 100644 --- a/examples/TUDELFT_V3_KITE/tutorial.py +++ b/examples/TUDELFT_V3_KITE/tutorial.py @@ -64,10 +64,11 @@ def main(): body_aero_CAD_CFD_polars_with_bridles = BodyAerodynamics.instantiate( n_panels=n_panels, # file_path=(cad_derived_geometry_dir / "config_kite_CAD_CFD_polars.yaml"), - file_path=(cad_derived_geometry_dir / "aero_geometry.yaml"), + file_path=(cad_derived_geometry_dir / "config_kite_CAD_CFD_polars.yaml"), spanwise_panel_distribution=spanwise_panel_distribution, - bridle_path=(cad_derived_geometry_dir / "struc_geometry.yaml"), - ml_models_dir=(Path(PROJECT_DIR) / "data" / "ml_models"), + bridle_path=( + cad_derived_geometry_dir / "struc_geometry_manually_adjusted.yaml" + ), ) body_aero_CAD_neuralfoil = BodyAerodynamics.instantiate( n_panels=n_panels, @@ -134,8 +135,6 @@ def main(): # / "interactive_plot.html", ) - breakpoint() - # Step 5: Plot polar curves for different angles of attack and side slip angles, and save results """ Compare the aerodynamic performance of different models by plotting lift, drag, and side force coefficients From 72439c8c9a6047d1b7f38566e6d56302dcf84aac Mon Sep 17 00:00:00 2001 From: ocayon Date: Tue, 28 Oct 2025 18:57:44 +0100 Subject: [PATCH 2/3] Change aoa 1/4 correction as optional bool in solver, false by default (#152) --- .../TUDELFT_V3_KITE/fit_polynomials_polars.py | 17 +- src/VSM/core/BodyAerodynamics.py | 7 +- src/VSM/core/Solver.py | 3 + src/VSM/plotting.py | 2 + src/VSM/stability_derivatives.py | 569 ++++++++++++++++++ 5 files changed, 585 insertions(+), 13 deletions(-) create mode 100644 src/VSM/stability_derivatives.py diff --git a/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py b/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py index 3d9157e..4577dd2 100644 --- a/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py +++ b/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py @@ -28,14 +28,14 @@ def main(): ### 1. defining paths PROJECT_DIR = Path(__file__).resolve().parents[2] - file_path = ( - Path(PROJECT_DIR) - / "data" - / "TUDELFT_V3_KITE" - / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" - ) - + # file_path = ( + # Path(PROJECT_DIR) + # / "data" + # / "TUDELFT_V3_KITE" + # / "CAD_derived_geometry" + # / "aero_geometry_CAD_CFD_polars.yaml" + # ) + file_path = Path(PROJECT_DIR) / "data" / "V9_KITE" / "config_kite.yaml" # bridle_path = ( # Path(PROJECT_DIR)s # / "data" @@ -50,7 +50,6 @@ def main(): solver = Solver(reference_point=[0, 0, 0]) ### 3. Loading kite geometry from CSV file and instantiating BodyAerodynamics - print(f"\nCreating corrected polar input with bridles") body_aero_polar_with_bridles = BodyAerodynamics.instantiate( n_panels=n_panels, file_path=file_path, diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index 92bfaec..e91290c 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -1024,6 +1024,7 @@ def compute_results( is_only_f_and_gamma_output, is_with_viscous_drag_correction, reference_point, + is_aoa_corrected, ): cl_array, cd_array, cm_array = ( @@ -1040,7 +1041,7 @@ def compute_results( drag = (cd_array * 0.5 * rho * Umag_array**2 * chord_array)[:, np.newaxis] moment = (cm_array * 0.5 * rho * Umag_array**2 * chord_array**2)[:, np.newaxis] - if aerodynamic_model_type == "VSM": + if is_aoa_corrected: alpha_corrected = self.update_effective_angle_of_attack_if_VSM( gamma_new, core_radius_fraction, @@ -1052,11 +1053,9 @@ def compute_results( ) alpha_uncorrected = alpha_array[:, np.newaxis] - elif aerodynamic_model_type == "LLT": + else: alpha_corrected = alpha_array[:, np.newaxis] alpha_uncorrected = alpha_array[:, np.newaxis] - else: - raise ValueError("Unknown aerodynamic model type, should be LLT or VSM") # Checking that va is not distributed input if len(self._va) != 3: raise ValueError("Calc.results not ready for va_distributed input") diff --git a/src/VSM/core/Solver.py b/src/VSM/core/Solver.py index bc3d85a..cb4a410 100644 --- a/src/VSM/core/Solver.py +++ b/src/VSM/core/Solver.py @@ -50,6 +50,7 @@ def __init__( artificial_damping: dict = {"k2": 0.1, "k4": 0.0}, is_with_simonet_artificial_viscosity: bool = False, simonet_artificial_viscosity_fva: float = None, + is_aoa_corrected: bool = False, ): """Initialize solver with configuration parameters. @@ -83,6 +84,7 @@ def __init__( self.is_only_f_and_gamma_output = is_only_f_and_gamma_output self.is_with_viscous_drag_correction = is_with_viscous_drag_correction self.reference_point = reference_point + self.is_aoa_corrected = is_aoa_corrected # === athmospheric properties === self.mu = mu self.rho = rho @@ -259,6 +261,7 @@ def solve(self, body_aero, gamma_distribution: np.ndarray = None) -> dict: self.is_only_f_and_gamma_output, self.is_with_viscous_drag_correction, self.reference_point, + self.is_aoa_corrected, ) return results diff --git a/src/VSM/plotting.py b/src/VSM/plotting.py index 2b939a1..423a78e 100644 --- a/src/VSM/plotting.py +++ b/src/VSM/plotting.py @@ -532,6 +532,8 @@ def plot_polars( polar_data[3] = df["CS"].values elif steering == "roll": polar_data[3] = np.degrees(np.arctan2(df["CS"], df["CL"])) + if angle_type == "angle_of_attack": + polar_data[3] = df["CL"].values / df["CD"].values if "CMx" in df.columns: polar_data[4] = df["CMx"].values if "CMy" in df.columns: diff --git a/src/VSM/stability_derivatives.py b/src/VSM/stability_derivatives.py new file mode 100644 index 0000000..a4f5dac --- /dev/null +++ b/src/VSM/stability_derivatives.py @@ -0,0 +1,569 @@ +"""Utilities for computing rigid-body aerodynamic stability derivatives. + +This module provides a helper function to evaluate force and moment +coefficient sensitivities with respect to kinematic angles (angle of attack, +sideslip) and body rotation rates (roll, pitch, yaw) by repeatedly invoking +the aerodynamic solver with finite-difference perturbations. +""" + +from __future__ import annotations + +from typing import Dict, Optional + +import numpy as np + + +CoeffDict = Dict[str, float] + + +def compute_rigid_body_stability_derivatives( + body_aero, + solver, + angle_of_attack: float, + side_slip: float, + velocity_magnitude: float, + roll_rate: float = 0.0, + pitch_rate: float = 0.0, + yaw_rate: float = 0.0, + step_sizes: Optional[Dict[str, float]] = None, + reference_point: Optional[np.ndarray] = None, + nondimensionalize_rates: bool = True, +) -> Dict[str, float]: + """Compute rigid-body stability derivatives for the current configuration. + + This function computes aerodynamic derivatives with respect to kinematic angles + (angle of attack, sideslip) and body rotation rates (roll, pitch, yaw). + + Parameters + ---------- + body_aero : VSM.core.BodyAerodynamics.BodyAerodynamics + Aerodynamic model instance that will be updated in-place. + solver : VSM.core.Solver.Solver + Solver configured for the analysis. + angle_of_attack : float + Baseline angle of attack in degrees. + side_slip : float + Baseline sideslip angle in degrees (positive for wind from the left/port side, + negative for wind from the right/starboard side). + velocity_magnitude : float + Magnitude of the freestream velocity (m/s). + roll_rate : float, optional + Baseline body roll rate ``p`` in rad/s. Defaults to 0.0. + pitch_rate : float, optional + Baseline body pitch rate ``q`` in rad/s. Defaults to 0.0. + yaw_rate : float, optional + Baseline body yaw rate ``r`` in rad/s. Defaults to 0.0. + step_sizes : dict, optional + Optional overrides for perturbation steps. Supported keys are + ``{"alpha", "beta", "p", "q", "r"}``. + Angle steps are in degrees (internally converted to radians for the derivative), + and rate steps are in rad/s. + reference_point : np.ndarray, optional + Reference point for moment calculation [x, y, z]. If None, defaults to + solver.reference_point if available, otherwise [0, 0, 0]. + nondimensionalize_rates : bool, optional + If True (default), rate derivatives are non-dimensionalized using: + hat_p = p * b / (2*V) + hat_q = q * c_MAC / (2*V) + hat_r = r * b / (2*V) + where b is wingspan, c_MAC is mean aerodynamic chord, and V is velocity magnitude. + This converts derivatives from per rad/s to per hat-rate (dimensionless). + If False, rate derivatives remain dimensional (per rad/s). + + Returns + ------- + dict + Dictionary with keys such as ``"dCx_dalpha"`` and ``"dCMz_dp"`` covering + stability derivatives with respect to angles and body rates: + + - Angle derivatives (per radian): dC*/dalpha, dC*/dbeta + - Rate derivatives: dC*/dp, dC*/dq, dC*/dr + - If nondimensionalize_rates=True: per hat-rate (dimensionless) + - If nondimensionalize_rates=False: per rad/s (dimensional) + + where C* ∈ {Cx, Cy, Cz, CMx, CMy, CMz} + + Notes + ----- + Derivatives are evaluated via central finite differences. Angular + sensitivities are returned per-radian. + + The reference_point parameter is critical for physically correct rotational + velocity calculations. The rotational velocity at any point r is computed as: + v_rot(r) = omega × (r - r_ref) + where omega is the body rate vector and r_ref is the reference point. + + Non-dimensionalization (when nondimensionalize_rates=True): + --------------------------------------------------------------- + Rate derivatives are converted to dimensionless form using: + - hat_p = p * b / (2*V) [roll rate] + - hat_q = q * c_MAC / (2*V) [pitch rate] + - hat_r = r * b / (2*V) [yaw rate] + + This converts derivatives from per rad/s to per hat-rate: + d(Coeff)/d(hat_p) = d(Coeff)/dp * (2*V/b) + d(Coeff)/d(hat_q) = d(Coeff)/dq * (2*V/c_MAC) + d(Coeff)/d(hat_r) = d(Coeff)/dr * (2*V/b) + + This is the standard convention used in flight dynamics and control texts, + making the derivatives directly compatible with 6-DOF equations of motion. + """ + + coeff_names = ("Cx", "Cy", "Cz", "CMx", "CMy", "CMz") + param_names = ("alpha", "beta", "p", "q", "r") + + default_steps = { + "alpha": np.rad2deg(0.005), # degrees + "beta": np.rad2deg(0.005), # degrees + "p": 0.01, # rad/s + "q": 0.01, # rad/s + "r": 0.01, # rad/s + } + if step_sizes: + for key, value in step_sizes.items(): + if key not in default_steps: + raise KeyError(f"Unsupported step key '{key}'. Allowed: {param_names}") + default_steps[key] = float(value) + + rates = {"p": roll_rate, "q": pitch_rate, "r": yaw_rate} + + # Use reference_point if provided, otherwise try solver.reference_point, else default to [0,0,0] + if reference_point is None: + if hasattr(solver, "reference_point"): + reference_point = np.array(solver.reference_point) + else: + reference_point = np.array([0.0, 0.0, 0.0]) + + def set_to_baseline() -> None: + body_aero.va_initialize( + Umag=velocity_magnitude, + angle_of_attack=angle_of_attack, + side_slip=side_slip, + yaw_rate=rates["r"], + pitch_rate=rates["q"], + roll_rate=rates["p"], + reference_point=reference_point, + ) + + def _solve_and_extract() -> CoeffDict: + results = solver.solve(body_aero) + va_vector = np.asarray(body_aero.va, dtype=float) + if va_vector.ndim != 1 or va_vector.size != 3: + raise ValueError("Expected a uniform apparent velocity vector of length 3.") + speed = np.linalg.norm(va_vector) + if speed <= 0.0: + raise ValueError( + "Freestream speed must be positive to compute derivatives." + ) + + q_inf = 0.5 * solver.rho * speed**2 + reference_area = results["projected_area"] + if reference_area <= 0.0: + raise ValueError("Reference area must be positive.") + + coeffs = { + "Cx": results["Fx"] / (q_inf * reference_area), + "Cy": results["Fy"] / (q_inf * reference_area), + "Cz": results["Fz"] / (q_inf * reference_area), + "CMx": results["cmx"], + "CMy": results["cmy"], + "CMz": results["cmz"], + } + return coeffs + + def _evaluate_with_angles(alpha_deg: float, beta_deg: float) -> CoeffDict: + body_aero.va_initialize( + Umag=velocity_magnitude, + angle_of_attack=alpha_deg, + side_slip=beta_deg, + yaw_rate=rates["r"], + pitch_rate=rates["q"], + roll_rate=rates["p"], + reference_point=reference_point, + ) + return _solve_and_extract() + + def _central_difference( + coeff_plus: CoeffDict, coeff_minus: CoeffDict, delta: float + ) -> CoeffDict: + if delta == 0.0: + raise ValueError("Delta for central difference cannot be zero.") + return { + name: (coeff_plus[name] - coeff_minus[name]) / (2.0 * delta) + for name in coeff_names + } + + derivatives: Dict[str, float] = {} + + # Baseline state + set_to_baseline() + + # Angle derivatives (per radian) + for param, step_key in (("alpha", "alpha"), ("beta", "beta")): + set_to_baseline() + step_deg = default_steps[step_key] + delta_rad = np.deg2rad(step_deg) + if param == "alpha": + coeff_plus = _evaluate_with_angles(angle_of_attack + step_deg, side_slip) + coeff_minus = _evaluate_with_angles(angle_of_attack - step_deg, side_slip) + else: + coeff_plus = _evaluate_with_angles(angle_of_attack, side_slip + step_deg) + coeff_minus = _evaluate_with_angles(angle_of_attack, side_slip - step_deg) + + diff = _central_difference(coeff_plus, coeff_minus, delta_rad) + + for coeff_name in coeff_names: + derivatives[f"d{coeff_name}_d{param}"] = diff[coeff_name] + + # Body rate derivatives (p, q, r) + for rate_name in ("p", "q", "r"): + set_to_baseline() + delta = default_steps[rate_name] + + # Create perturbed rate dictionaries + rates_plus = rates.copy() + rates_minus = rates.copy() + rates_plus[rate_name] += delta + rates_minus[rate_name] -= delta + + # Evaluate with positive perturbation + body_aero.va_initialize( + Umag=velocity_magnitude, + angle_of_attack=angle_of_attack, + side_slip=side_slip, + yaw_rate=rates_plus["r"], + pitch_rate=rates_plus["q"], + roll_rate=rates_plus["p"], + reference_point=reference_point, + ) + coeff_plus = _solve_and_extract() + + # Evaluate with negative perturbation + body_aero.va_initialize( + Umag=velocity_magnitude, + angle_of_attack=angle_of_attack, + side_slip=side_slip, + yaw_rate=rates_minus["r"], + pitch_rate=rates_minus["q"], + roll_rate=rates_minus["p"], + reference_point=reference_point, + ) + coeff_minus = _solve_and_extract() + + diff = _central_difference(coeff_plus, coeff_minus, delta) + + for coeff_name in coeff_names: + derivatives[f"d{coeff_name}_d{rate_name}"] = diff[coeff_name] + + # Restore baseline state before returning + set_to_baseline() + + # Non-dimensionalize rate derivatives if requested + if nondimensionalize_rates: + # Get geometric properties from baseline solution + results_baseline = solver.solve(body_aero) + b = results_baseline["wing_span"] # wingspan + + # Compute mean aerodynamic chord from projected area and span + # Note: c_MAC = S / b is a reasonable approximation for simple geometries + # For more complex planforms, a weighted average chord should be computed + S = results_baseline["projected_area"] + c_mac = S / b + print( + f"\n --> Computed c_MAC = {c_mac:.3f} m from S={S:.3f} m² and b={b:.3f} m." + ) + + V = velocity_magnitude + + # Scaling factors to convert from per rad/s to per hat-rate: + # hat_p = p * b / (2*V) --> d()/dp [per rad/s] * (2*V/b) = d()/d(hat_p) + # hat_q = q * c_MAC / (2*V) --> d()/dq [per rad/s] * (2*V/c_MAC) = d()/d(hat_q) + # hat_r = r * b / (2*V) --> d()/dr [per rad/s] * (2*V/b) = d()/d(hat_r) + scale_p = (2.0 * V) / b + scale_q = (2.0 * V) / c_mac + scale_r = (2.0 * V) / b + + # Apply scaling to all rate derivatives + for coeff_name in coeff_names: + key_p = f"d{coeff_name}_dp" + key_q = f"d{coeff_name}_dq" + key_r = f"d{coeff_name}_dr" + + if key_p in derivatives: + derivatives[key_p] *= scale_p # now per hat_p + if key_q in derivatives: + derivatives[key_q] *= scale_q # now per hat_q + if key_r in derivatives: + derivatives[key_r] *= scale_r # now per hat_r + + return derivatives + + +def map_derivatives_to_aircraft_frame( + derivatives_vsm: Dict[str, float], +) -> Dict[str, float]: + """ + Transform stability derivatives from VSM frame (x rearward, y right, z up) + to aircraft frame (x forward, y right, z down) using the proper rotation + R_y(pi) = diag(-1, 1, -1). + + Signs: + s(Cx)=s(Cz)=s(CMx)=s(CMz) = -1 + s(Cy)=s(CMy) = +1 + + Mapping of independent variables: + alpha' = alpha + beta' = -beta + p' = -p, q' = q, r' = -r + """ + print("\n" + "=" * 60) + print("REFERENCE FRAME TRANSFORMATION") + print("=" * 60) + + print("\nVSM Reference Frame (used above):") + print(" x: rearward (LE → TE)") + print(" y: right wing") + print(" z: upward") + print(" β: positive for wind from LEFT (port, +Vy)") + print(" α: positive for nose up") + print(" Body rates (right-hand rule with z-up):") + print(" p (roll): positive = LEFT wing DOWN (CCW looking forward)") + print(" q (pitch): positive = nose UP (CW looking from right wing)") + print(" r (yaw): positive = nose LEFT (CCW looking down)") + + print("\n" + "-" * 40) + print(f"Aircraft Reference Frame (standard aerospace): ") + print(f" x: forward (tail → nose)") + print(f" y: right wing") + print(f" z: downward") + print(f" Mapping used: R_y(pi) = diag(-1, 1, -1) [proper 180° about +y]") + print(f" Implications:") + print(f" - Forces: (Cx', Cy', Cz') = (-Cx, Cy, -Cz)") + print(f" - Moments: (CMx',CMy',CMz')= (-CMx, CMy, -CMz)") + print(f" - Angles: alpha' = alpha, beta' = -beta") + print(f" - Rates: (p', q', r') = (-p, q, -r)") + + print("\n" + "=" * 60) + print("STABILITY DERIVATIVES IN AIRCRAFT FRAME") + print("=" * 60) + s = {"Cx": -1, "Cy": +1, "Cz": -1, "CMx": -1, "CMy": +1, "CMz": -1} + + out: Dict[str, float] = {} + + # Angle derivatives + for coeff, sc in s.items(): + ka = f"d{coeff}_dalpha" # alpha' = alpha (no flip) + kb = f"d{coeff}_dbeta" # beta' = -beta (extra minus) + + if ka in derivatives_vsm: + out[ka] = sc * derivatives_vsm[ka] + if kb in derivatives_vsm: + out[kb] = -sc * derivatives_vsm[kb] + + # Rate derivatives + for coeff, sc in s.items(): + kp = f"d{coeff}_dp" # p' = -p (extra minus) + kq = f"d{coeff}_dq" # q' = q (no flip) + kr = f"d{coeff}_dr" # r' = -r (extra minus) + + if kp in derivatives_vsm: + out[kp] = -sc * derivatives_vsm[kp] + if kq in derivatives_vsm: + out[kq] = sc * derivatives_vsm[kq] + if kr in derivatives_vsm: + out[kr] = -sc * derivatives_vsm[kr] + + return out + + +import numpy as np +import math +from typing import Dict, Tuple, Optional, Sequence + + +# --- Reuse the adapter that maps aircraft inputs -> VSM solver and back ------------- +def compute_aircraft_frame_coeffs_from_solver( + body_aero, + solver, + alpha_rad: float, + beta_rad: float, + V: float, + p: float = 0.0, + q: float = 0.0, + r: float = 0.0, + reference_point: Optional[np.ndarray] = None, +) -> Tuple[Dict[str, float], float, float]: + # map aircraft->VSM + alpha_deg_vsm = math.degrees(alpha_rad) + beta_deg_vsm = -math.degrees(beta_rad) + p_vsm, q_vsm, r_vsm = -p, q, -r + + if reference_point is None: + if hasattr(solver, "reference_point"): + reference_point = np.array(solver.reference_point, dtype=float) + else: + reference_point = np.array([0.0, 0.0, 0.0], dtype=float) + + body_aero.va_initialize( + Umag=V, + angle_of_attack=alpha_deg_vsm, + side_slip=beta_deg_vsm, + yaw_rate=r_vsm, + pitch_rate=q_vsm, + roll_rate=p_vsm, + reference_point=reference_point, + ) + results = solver.solve(body_aero) + + va_vec = np.asarray(body_aero.va, float) + Va = float(np.linalg.norm(va_vec)) + if Va <= 0: + raise ValueError("Speed must be >0") + rho = float(getattr(solver, "rho", 1.225)) + q_inf = 0.5 * rho * Va * Va + S = float(results["projected_area"]) + if S <= 0: + raise ValueError("projected_area must be >0") + + Fx, Fy, Fz = float(results["Fx"]), float(results["Fy"]), float(results["Fz"]) + CMx, CMy, CMz = float(results["cmx"]), float(results["cmy"]), float(results["cmz"]) + + # VSM -> aircraft mapping: R_y(pi) = diag(-1, +1, -1) + Cx_vsm = Fx / (q_inf * S) + Cy_vsm = Fy / (q_inf * S) + Cz_vsm = Fz / (q_inf * S) + coeffs_air = { + "CX": -Cx_vsm, + "CY": Cy_vsm, + "CZ": -Cz_vsm, + "Cl": -CMx, + "Cm": CMy, + "Cn": -CMz, + } + return coeffs_air, q_inf, S + + +# --- Utilities --------------------------------------------------------------------- +def fit_quad_alpha( + alpha_vec: np.ndarray, y_alpha: np.ndarray +) -> Tuple[float, float, float]: + """Least-squares fit y(α) ≈ c2 α^2 + c1 α + c0; α in radians.""" + A = np.column_stack([alpha_vec**2, alpha_vec, np.ones_like(alpha_vec)]) + c2, c1, c0 = np.linalg.lstsq(A, y_alpha, rcond=None)[0] + return float(c2), float(c1), float(c0) + + +def central_diff(f_plus: float, f_minus: float, h: float) -> float: + return (f_plus - f_minus) / (2.0 * h) + + +# --- Main builder ------------------------------------------------------------------ +def build_malz_coeff_table_from_solver( + body_aero, + solver, + V: float, + b: float, + c: float, + alpha_grid_deg: Sequence[float] = tuple(np.linspace(-15, 15, 13)), + beta_step_deg: float = 1.0, + p_hat_step: float = 0.02, # small dimensionless hat-rate steps + q_hat_step: float = 0.02, + r_hat_step: float = 0.02, + reference_point: Optional[np.ndarray] = None, +) -> Dict[str, Tuple[float, float, float]]: + """ + Returns a dict like COEF with keys: + CX0,CXb,CXp,CXq,CXr, CY0,CYb,CYp,CYq,CYr, CZ0,CZb,CZp,CZq,CZr, + Cl0,Clb,Clp,Clq,Clr, Cm0,Cmb,Cmp,Cmq,Cmr, Cn0,Cnb,Cnp,Cnq,Cnr + Each value is a (c2,c1,c0) tuple per Eq. (19), α in radians. + + All evaluations are at controls = 0. + """ + alpha_grid = np.radians(np.array(alpha_grid_deg, float)) + # small angle/rate steps + dbeta = math.radians(beta_step_deg) + + # convert hat-rate steps to actual rad/s perturbations at this V + # p̂ = b p / (2 V) => p = (2 V / b) p̂ ; etc. + p_step = (2.0 * V / b) * p_hat_step + q_step = (2.0 * V / c) * q_hat_step + r_step = (2.0 * V / b) * r_hat_step + + # storage per α for each output we will fit + names = ["CX", "CY", "CZ", "Cl", "Cm", "Cn"] + slices_0 = {k: [] for k in names} + slices_db = {k: [] for k in names} + slices_dp = {k: [] for k in names} + slices_dq = {k: [] for k in names} + slices_dr = {k: [] for k in names} + + beta0 = 0.0 + p0 = q0 = r0 = 0.0 + + for a in alpha_grid: + # baseline + c0, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0, q0, r0, reference_point + ) + for k in names: + slices_0[k].append(c0[k]) + + # beta derivative at this α + cp, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0 + dbeta, V, p0, q0, r0, reference_point + ) + cm, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0 - dbeta, V, p0, q0, r0, reference_point + ) + for k in names: + slices_db[k].append(central_diff(cp[k], cm[k], dbeta)) + + # p̂ derivative (via p step in rad/s) + cp_, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0 + p_step, q0, r0, reference_point + ) + cm_, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0 - p_step, q0, r0, reference_point + ) + for k in names: + # convert from per (rad/s) to per hat_p by multiplying (2V/b) + slices_dp[k].append(central_diff(cp_[k], cm_[k], p_step) * (2.0 * V / b)) + + # q̂ derivative + cq_, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0, q0 + q_step, r0, reference_point + ) + cmq, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0, q0 - q_step, r0, reference_point + ) + for k in names: + slices_dq[k].append(central_diff(cq_[k], cmq[k], q_step) * (2.0 * V / c)) + + # r̂ derivative + cr_, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0, q0, r0 + r_step, reference_point + ) + cmr, _, _ = compute_aircraft_frame_coeffs_from_solver( + body_aero, solver, a, beta0, V, p0, q0, r0 - r_step, reference_point + ) + for k in names: + slices_dr[k].append(central_diff(cr_[k], cmr[k], r_step) * (2.0 * V / b)) + + # Fit each α-slice to quadratic c2 α^2 + c1 α + c0, α in radians + COEF = {} + + # Baseline (*0) + for k in names: + COEF[k + "0"] = fit_quad_alpha(alpha_grid, np.array(slices_0[k])) + + # Beta (*b) + for k in names: + COEF[k + "b"] = fit_quad_alpha(alpha_grid, np.array(slices_db[k])) + + # Rates (*p,*q,*r) + for k in names: + COEF[k + "p"] = fit_quad_alpha(alpha_grid, np.array(slices_dp[k])) + COEF[k + "q"] = fit_quad_alpha(alpha_grid, np.array(slices_dq[k])) + COEF[k + "r"] = fit_quad_alpha(alpha_grid, np.array(slices_dr[k])) + + return COEF From eae1cb6b9776e9a1c0bb212d93ad81c0e7b64e42 Mon Sep 17 00:00:00 2001 From: ocayon Date: Tue, 28 Oct 2025 19:07:27 +0100 Subject: [PATCH 3/3] #152 add test of correction --- .../TUDELFT_V3_KITE/aoa_correction_test.py | 140 ++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 examples/TUDELFT_V3_KITE/aoa_correction_test.py diff --git a/examples/TUDELFT_V3_KITE/aoa_correction_test.py b/examples/TUDELFT_V3_KITE/aoa_correction_test.py new file mode 100644 index 0000000..fb2ede2 --- /dev/null +++ b/examples/TUDELFT_V3_KITE/aoa_correction_test.py @@ -0,0 +1,140 @@ +from pathlib import Path +import numpy as np +from VSM.core.BodyAerodynamics import BodyAerodynamics +from VSM.core.Solver import Solver +from VSM.plotting import ( + plot_polars, +) +from VSM.plot_geometry_matplotlib import plot_geometry +from VSM.plot_geometry_plotly import interactive_plot + + +def main(): + """ + Example: 3D Aerodynamic Analysis of TUDELFT_V3_KITE using VSM + + This script demonstrates the workflow for performing a 3D aerodynamic analysis of the TUDELFT_V3_KITE + using the Vortex Step Method (VSM) library. The workflow is structured as follows: + + Step 1: Instantiate BodyAerodynamics objects from different YAML configuration files. + - Each YAML config defines the geometry and airfoil/polar data for a specific modeling approach. + - Supported approaches include: + - Breukels regression (empirical model) + - CFD-based polars + - NeuralFoil-based polars + - Masure regression (machine learning model) + - Inviscid theory + + Step 2: Set inflow conditions for each aerodynamic object. + - Specify wind speed (Umag), angle of attack, side slip, and yaw rate. + - Initialize the apparent wind for each BodyAerodynamics object. + + Step 3: Plot the kite geometry using Matplotlib. + - Visualize the panel mesh, control points, and aerodynamic centers. + + Step 4: Create an interactive 3D plot using Plotly. + - Allows for interactive exploration of the geometry and panel arrangement. + + Step 5: Plot and save polar curves for different angles of attack and side slip angles. + - Compare the results of different aerodynamic models. + - Optionally include literature/CFD data for validation. + + Step 5a: Plot alpha sweep (angle of attack variation). + Step 5b: Plot beta sweep (side slip variation). + + Returns: + None + """ + ### 1. defining paths + PROJECT_DIR = Path(__file__).resolve().parents[2] + + ### 2. defining settings + n_panels = 50 + spanwise_panel_distribution = "uniform" + solver_base_version = Solver(reference_point=np.array([0.0, 0.0, 0.0])) + solver_aoa_correction = Solver( + reference_point=np.array([0.0, 0.0, 0.0]), is_aoa_corrected=True + ) + + # Step 1: Instantiate BodyAerodynamics objects from different YAML configs + cad_derived_geometry_dir = ( + Path(PROJECT_DIR) / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" + ) + body_aero_masure_regression = BodyAerodynamics.instantiate( + n_panels=n_panels, + file_path=( + cad_derived_geometry_dir / "aero_geometry_CAD_masure_regression.yaml" + ), + ml_models_dir=(Path(PROJECT_DIR) / "data" / "ml_models"), + spanwise_panel_distribution=spanwise_panel_distribution, + ) + + # Step 2: Set inflow conditions for each aerodynamic object + """ + Set the wind speed, angle of attack, side slip, and yaw rate for each BodyAerodynamics object. + This initializes the apparent wind vector and prepares the objects for analysis. + """ + Umag = 3.15 + angle_of_attack = 6.8 + side_slip = 0 + yaw_rate = 0 + + body_aero_masure_regression.va_initialize( + Umag, angle_of_attack, side_slip, yaw_rate + ) + + # Step 5: Plot polar curves for different angles of attack and side slip angles, and save results + """ + Compare the aerodynamic performance of different models by plotting lift, drag, and side force coefficients + as a function of angle of attack (alpha sweep) and side slip (beta sweep). + Literature/CFD data can be included for validation. + """ + save_folder = Path(PROJECT_DIR) / "results" / "TUDELFT_V3_KITE" + + # Step 5a: Plot alpha sweep (angle of attack) + path_cfd_lebesque_alpha = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "CFD_RANS_Rey_10e5_Poland2025_alpha_sweep_beta_0.csv" + ) + path_wt_data = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "V3_CL_CD_CS_alpha_sweep_for_beta_0_WindTunnel_Poland_2025_Rey_560e4.csv" + ) + plot_polars( + solver_list=[ + solver_aoa_correction, + solver_base_version, + ], + body_aero_list=[ + body_aero_masure_regression, + body_aero_masure_regression, + ], + label_list=[ + "VSM Aoa Correction", + "VSM No Aoa Correction", + "CFD_RANS_Rey_10e5_Poland2025_alpha_sweep_beta_0", + "Wind Tunnel Data Poland 2025", + ], + literature_path_list=[path_cfd_lebesque_alpha, path_wt_data], + angle_range=np.linspace(0, 20, 21), + angle_type="angle_of_attack", + angle_of_attack=0, + side_slip=0, + yaw_rate=0, + Umag=Umag, + title="alphasweep", + data_type=".pdf", + save_path=Path(save_folder), + is_save=True, + is_show=False, + ) + + +if __name__ == "__main__": + main()