diff --git a/README.md b/README.md index 6456dd0..d2178d8 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,34 @@ # Vortex Step Method -The Vortex Step Method (VSM) is an enhanced lifting line method that improves upon the classic approach by solving the circulation system at the three-quarter chord position. This adjustment allows for more accurate calculations of lift and drag forces, particularly addressing the shortcomings in induced drag prediction. VSM is further refined by coupling it with 2D viscous airfoil polars, making it well-suited for complex geometries, including low aspect ratio wings, as well as configurations with sweep, dihedral, and anhedral angles, typical for leading-edge inflatable (LEI) kites, that are used in airborne wind energy production, boat-towing and kite-surfing. An open-source example kite is the [TU Delft V3 Kite](https://awegroup.github.io/TUDELFT_V3_KITE/), of which a video is shown below, made using the internal plotting library. +The Vortex Step Method (VSM) is an enhanced lifting line method that improves upon the classic approach by solving the circulation system at the three-quarter chord position. This adjustment allows for more accurate calculations of lift and drag forces, particularly addressing the shortcomings in induced drag prediction. VSM is further refined by coupling it with 2D viscous airfoil polars, making it well-suited for complex geometries, including low aspect ratio wings, as well as configurations with sweep, dihedral, and anhedral angles, typical for leading-edge inflatable (LEI) kites, that are used in airborne wind energy production, boat-towing and kite-surfing. -![](docs/TUDELFT_V3_KITE_plotly.gif) +A Julia version of this project is available at [VortexStepMethod.jl](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) -The software presented here includes examples for: a rectangular wing and a leading-edge inflatable kite. -A Julia version of this project is available at [VortexStepMethod.jl](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl) +## Reference Frame + +The VSM uses a body-fixed reference frame with the following conventions: + +**Coordinate Axes:** +- **x-axis**: Rearward (leading edge → trailing edge) +- **y-axis**: Right wing (when looking from behind the kite) +- **z-axis**: Upward (from wing tip to mid-span) + +**Aerodynamic Angles:** +- **Angle of Attack (α)**: Positive for nose up +- **Sideslip (β)**: Positive for wind from the left/port side (positive Vy) + +**Body Rotation Rates** (right-hand rule): +- **Roll rate (p)**: Positive for left wing down (counter-clockwise looking forward) +- **Pitch rate (q)**: Positive for nose up (clockwise looking from right wing) +- **Yaw rate (r)**: Positive for nose left (counter-clockwise looking down) + +The reference frame is illustrated below for the open-source example kite, the [TU Delft V3 Kite](https://awegroup.github.io/TUDELFT_V3_KITE/). + + +![](docs/TUDELFT_V3_KITE_reference_frame.png) + +**Aircraft Frame Transformation:** For stability derivatives in standard aircraft coordinates (x-forward, y-right, z-down), use `map_derivatives_to_aircraft_frame()` from the `VSM.stability_derivatives` module. See `examples/TUDELFT_V3_KITE/evaluate_stability_derivatives.py` for usage. + ## Key Features @@ -99,41 +122,6 @@ For detailed documentation, please refer to the following resources. deactivate ``` -## Quick Start - -Here's a minimal example to get started: - -```python -from pathlib import Path -from VSM.core.BodyAerodynamics import BodyAerodynamics -from VSM.core.Solver import Solver - -# Load kite geometry from YAML configuration -config_path = "data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_CFD_polars.yaml" -body_aero = BodyAerodynamics.instantiate( - n_panels=30, - file_path=config_path, - spanwise_panel_distribution="uniform" -) - -# Set flow conditions -body_aero.va_initialize( - Umag=10.0, # Velocity magnitude [m/s] - angle_of_attack=6.0, # Angle of attack [deg] - side_slip=0.0, # Sideslip angle [deg] -) - -# Solve and get results -solver = Solver() -results = solver.solve(body_aero) - -print(f"CL = {results['cl']:.3f}") -print(f"CD = {results['cd']:.3f}") -print(f"L/D = {results['cl']/results['cd']:.2f}") -``` - -For more examples, see the `examples/` directory. - ## Dependencies - numpy - Numerical computing @@ -182,9 +170,7 @@ See individual files for detailed documentation and usage instructions. For large-scale parametric studies: 1. **Use fewer panels during exploration** (`n_panels=20-30`), then increase for final results -2. **Numba JIT compilation** automatically accelerates vortex calculations (first run slower) -3. **Parallel studies**: Run multiple configurations using `multiprocessing` or `joblib` -4. **Cache geometries**: Instantiate `BodyAerodynamics` once, reuse with different flow conditions +2. **Cache geometries**: Instantiate `BodyAerodynamics` once, reuse with different flow conditions Example: ```python diff --git a/changelog.md b/changelog.md index 8a500b7..86ab928 100644 --- a/changelog.md +++ b/changelog.md @@ -5,6 +5,15 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). +## [2.0.1] - 21-10-2025 + +### Added +- `evaluate_stability_derivatives.py` example script demonstrating the use of the new + `compute_rigid_body_stability_derivatives()` function for calculating aerodynamic stability derivatives, prints results at trim conditions in the kite reference frame (x-backward LE to TE, y-right from kite's perspective, z-up) and in the aircraft convention (x-forward, y-right, z-down). + +### Fixed +- Corrected a misplaced minus sign in the sideslip handling inside `BodyAerodynamics.va_initialize`, restoring the intended sign convention for moment predictions and β-related stability derivatives. + ## [2.0.0] - 08-10-2025 ### ⚠️ Breaking Changes diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_CFD_polars.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_CFD_polars.yaml similarity index 100% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_CFD_polars.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_CFD_polars.yaml diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_airfoil.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_airfoil.yaml similarity index 70% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_airfoil.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_airfoil.yaml index 230c0ab..fe37a84 100644 --- a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_airfoil.yaml +++ b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_airfoil.yaml @@ -101,129 +101,3 @@ wing_airfoils: - [17, neuralfoil, {dat_file_path: "clark_y.dat", model_size: "xxxlarge", xtr_lower: 0.0, xtr_upper: 1.0, n_crit: 9.0}] - [18, neuralfoil, {dat_file_path: "clark_y.dat", model_size: "xxxlarge", xtr_lower: 0.0, xtr_upper: 1.0, n_crit: 9.0}] - [19, neuralfoil, {dat_file_path: "clark_y.dat", model_size: "xxxlarge", xtr_lower: 0.0, xtr_upper: 1.0, n_crit: 9.0}] - - - -bridle_nodes: - # --------------------------------------------------------------- - # headers: - # - id: integer, unique identifier for the node - # - x: x-coordinate [m] - # - y: y-coordinate [m] - # - z: z-coordinate [m] - # - type: node type, either 'knot' or 'pulley' - # --------------------------------------------------------------- - headers: [id, x, y, z, type] - data: - - [21, -0.67653494, 2.10580567, 7.08413599, knot] - - [22, -0.79196791, 1.01882151, 7.57069579, knot] - - [23, -0.79196791, -1.01882151, 7.57069579, knot] - - [24, -0.67653494, -2.10580567, 7.08413599, knot] - - [25, -0.46444477, 0.77505128, 4.23316161, knot] - - [26, -0.46444477, -0.77505128, 4.23316161, knot] - - [27, 0.54344694, 1.57655435, 5.91086207, knot] - - [28, 0.62960591, 0.92728966, 6.64639402, knot] - - [29, 0.62960591, -0.92728966, 6.64639402, knot] - - [30, 0.54344694, -1.57655435, 5.91086207, knot] - - [31, 0.39092886, 0.71221536, 4.26118728, knot] - - [32, 0.39092886, -0.71221536, 4.26118728, knot] - - [33, 0.30679124, 0.49235914, 3.39885672, pulley] - - [34, 0.28625931, 0.0, 3.19030736, knot] - - [35, 0.30679124, -0.49235914, 3.39885672, pulley] - - [36, 0.13269336, 0.23394851, 1.48060245, knot] - - [37, 0.13269336, -0.23394851, 1.48060245, knot] - -bridle_lines: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - rest_length: measured rest length [m] - # - diameter: line diameter [m] - # - material: string, material type (e.g., dyneema) - # - density: material density [kg/m^3] - # --------------------------------------------------------------- - headers: [name, rest_length, diameter, material, density] - data: - # front lines - - [a6, 11.790, 0.002, dyneema,970] - - [A5, 2.870, 0.002, dyneema,970] - - [A3, 3.450, 0.002, dyneema,970] - - [A2, 3.500, 0.002, dyneema,970] - - [A1, 3.690, 0.002, dyneema,970] - - [AII, 3.250, 0.002, dyneema,970] - - [AI, 3.620, 0.002, dyneema,970] - - [amain, 3.910, 0.002, dyneema,970] - # rear lines - - [br1, 4.410, 0.002, dyneema,970] - - [br2, 4.170, 0.002, dyneema,970] - - [br3, 4.140, 0.002, dyneema,970] - - [br5, 3.580, 0.002, dyneema,970] - - [br6, 13.560, 0.002, dyneema,970] - - [BRI, 2.360, 0.002, dyneema,970] - - [BRII, 2.080, 0.002, dyneema,970] - - [BR-main-1, .700, 0.002, dyneema,970] - - [M-line, 2.410, 0.002, dyneema,970] - # - [S, 1.730, 0.002, dyneema,970] # measured - # - [Power Tape, 3.330, 0.002, polyester] # measured - - [Steering Tape, 1.530, 0.002, dyneema,970] # trying 0.2m smalller - - [Power Tape, 2.900, 0.002, dyneema,970] # trying 0.2m smaller - -bridle_connections: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - ci: integer, node id (start) - # - cj: integer, node id (end) - # - ck: integer, third node id (only for pulleys, else omitted or 0) - # --------------------------------------------------------------- - headers: [name, ci, cj, ck] - data: - # front lines - - [a6,1,25,33] #e.g. this is a pulley - - [A5,3,21] - - [A3,5,21] - - [A2,7,22] - - [A1,9,22] - - [A1,11,23] - - [A2,13,23] - - [A3,15,24] - - [A5,17,24] - - [a6,19,26,35] - - [AII,21,25] - - [AI,22,25] - - [AI,23,26] - - [AII,24,26] - - [amain,25,0] - - [amain,26,0] - # rear lines - - [a6,1,33,25] - - [br6,2,33,36] - - [br6,2,36,33] - - [br5,4,27] - - [br3,6,27] - - [br2,8,28] - - [br1,10,28] - - [br1,12,29] - - [br2,14,29] - - [br3,16,30] - - [br5,18,30] - - [a6,19,35,26] - - [br6,20,35,37] - - [br6,20,37,35] - - [BRII,27,31] - - [BRI,28,31] - - [BRI,29,32] - - [BRII,30,32] - - [BR-main-1,31,33] - - [BR-main-1,32,35] - - [M-line,33,36,34] - - [M-line,33,34,36] - - [M-line,35,34,37] - - [M-line,35,37,34] - - [Steering Tape,36,0] - - [Power Tape,34,0] - - [Steering Tape,37,0] - - - - diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_breukels_regression.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_breukels_regression.yaml similarity index 69% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_breukels_regression.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_breukels_regression.yaml index 4fb2207..59b463b 100644 --- a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_breukels_regression.yaml +++ b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_breukels_regression.yaml @@ -104,123 +104,3 @@ wing_airfoils: - [19,breukels_regression, {t: 0.07, eta: 0.02, kappa: 0.59, delta: -0.2, lambda: 0.29, phi: 0.07}] - -bridle_nodes: - # --------------------------------------------------------------- - # headers: - # - id: integer, unique identifier for the node - # - x: x-coordinate [m] - # - y: y-coordinate [m] - # - z: z-coordinate [m] - # - type: node type, either 'knot' or 'pulley' - # --------------------------------------------------------------- - headers: [id, x, y, z, type] - data: - - [21, -0.67653494, 2.10580567, 7.08413599, knot] - - [22, -0.79196791, 1.01882151, 7.57069579, knot] - - [23, -0.79196791, -1.01882151, 7.57069579, knot] - - [24, -0.67653494, -2.10580567, 7.08413599, knot] - - [25, -0.46444477, 0.77505128, 4.23316161, knot] - - [26, -0.46444477, -0.77505128, 4.23316161, knot] - - [27, 0.54344694, 1.57655435, 5.91086207, knot] - - [28, 0.62960591, 0.92728966, 6.64639402, knot] - - [29, 0.62960591, -0.92728966, 6.64639402, knot] - - [30, 0.54344694, -1.57655435, 5.91086207, knot] - - [31, 0.39092886, 0.71221536, 4.26118728, knot] - - [32, 0.39092886, -0.71221536, 4.26118728, knot] - - [33, 0.30679124, 0.49235914, 3.39885672, pulley] - - [34, 0.28625931, 0.0, 3.19030736, knot] - - [35, 0.30679124, -0.49235914, 3.39885672, pulley] - - [36, 0.13269336, 0.23394851, 1.48060245, knot] - - [37, 0.13269336, -0.23394851, 1.48060245, knot] - -bridle_lines: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - rest_length: measured rest length [m] - # - diameter: line diameter [m] - # - material: string, material type (e.g., dyneema) - # - density: material density [kg/m^3] - # --------------------------------------------------------------- - headers: [name, rest_length, diameter, material, density] - data: - # front lines - - [a6, 11.790, 0.002, dyneema,970] - - [A5, 2.870, 0.002, dyneema,970] - - [A3, 3.450, 0.002, dyneema,970] - - [A2, 3.500, 0.002, dyneema,970] - - [A1, 3.690, 0.002, dyneema,970] - - [AII, 3.250, 0.002, dyneema,970] - - [AI, 3.620, 0.002, dyneema,970] - - [amain, 3.910, 0.002, dyneema,970] - # rear lines - - [br1, 4.410, 0.002, dyneema,970] - - [br2, 4.170, 0.002, dyneema,970] - - [br3, 4.140, 0.002, dyneema,970] - - [br5, 3.580, 0.002, dyneema,970] - - [br6, 13.560, 0.002, dyneema,970] - - [BRI, 2.360, 0.002, dyneema,970] - - [BRII, 2.080, 0.002, dyneema,970] - - [BR-main-1, .700, 0.002, dyneema,970] - - [M-line, 2.410, 0.002, dyneema,970] - # - [S, 1.730, 0.002, dyneema,970] # measured - # - [Power Tape, 3.330, 0.002, polyester] # measured - - [Steering Tape, 1.530, 0.002, dyneema,970] # trying 0.2m smalller - - [Power Tape, 2.900, 0.002, dyneema,970] # trying 0.2m smaller - -bridle_connections: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - ci: integer, node id (start) - # - cj: integer, node id (end) - # - ck: integer, third node id (only for pulleys, else omitted or 0) - # --------------------------------------------------------------- - headers: [name, ci, cj, ck] - data: - # front lines - - [a6,1,25,33] #e.g. this is a pulley - - [A5,3,21] - - [A3,5,21] - - [A2,7,22] - - [A1,9,22] - - [A1,11,23] - - [A2,13,23] - - [A3,15,24] - - [A5,17,24] - - [a6,19,26,35] - - [AII,21,25] - - [AI,22,25] - - [AI,23,26] - - [AII,24,26] - - [amain,25,0] - - [amain,26,0] - # rear lines - - [a6,1,33,25] - - [br6,2,33,36] - - [br6,2,36,33] - - [br5,4,27] - - [br3,6,27] - - [br2,8,28] - - [br1,10,28] - - [br1,12,29] - - [br2,14,29] - - [br3,16,30] - - [br5,18,30] - - [a6,19,35,26] - - [br6,20,35,37] - - [br6,20,37,35] - - [BRII,27,31] - - [BRI,28,31] - - [BRI,29,32] - - [BRII,30,32] - - [BR-main-1,31,33] - - [BR-main-1,32,35] - - [M-line,33,36,34] - - [M-line,33,34,36] - - [M-line,35,34,37] - - [M-line,35,37,34] - - [Steering Tape,36,0] - - [Power Tape,34,0] - - [Steering Tape,37,0] diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_masure_regression.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_masure_regression.yaml similarity index 54% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_masure_regression.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_masure_regression.yaml index 3f15f1b..83a1de4 100644 --- a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_masure_regression.yaml +++ b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_masure_regression.yaml @@ -83,144 +83,24 @@ wing_airfoils: reynolds: !!float 1e6 # Reynolds number headers: [airfoil_id, type, info_dict] data: - - [1, masure_regression, {t: 0.07, eta: 0.2, kappa: 0.95, delta: -2, lambda: 0.65, phi: 0.25}] - - [2, masure_regression, {t: 0.07, eta: 0.19, kappa: 0.93, delta: -1.9, lambda: 0.63, phi: 0.24}] - - [3, masure_regression, {t: 0.07, eta: 0.18, kappa: 0.91, delta: -1.8, lambda: 0.61, phi: 0.23}] - - [4, masure_regression, {t: 0.07, eta: 0.17, kappa: 0.89, delta: -1.7, lambda: 0.59, phi: 0.22}] - - [5, masure_regression, {t: 0.07, eta: 0.16, kappa: 0.87, delta: -1.6, lambda: 0.57, phi: 0.21}] - - [6, masure_regression, {t: 0.07, eta: 0.15, kappa: 0.85, delta: -1.5, lambda: 0.55, phi: 0.20}] - - [7, masure_regression, {t: 0.07, eta: 0.14, kappa: 0.83, delta: -1.4, lambda: 0.53, phi: 0.19}] - - [8, masure_regression, {t: 0.07, eta: 0.13, kappa: 0.81, delta: -1.3, lambda: 0.51, phi: 0.18}] - - [9, masure_regression, {t: 0.07, eta: 0.12, kappa: 0.79, delta: -1.2, lambda: 0.49, phi: 0.17}] - - [10, masure_regression, {t: 0.07, eta: 0.11, kappa: 0.77, delta: -1.1, lambda: 0.47, phi: 0.16}] - - [11, masure_regression, {t: 0.07, eta: 0.10, kappa: 0.75, delta: -1.0, lambda: 0.45, phi: 0.15}] - - [12, masure_regression, {t: 0.07, eta: 0.09, kappa: 0.73, delta: -0.9, lambda: 0.43, phi: 0.14}] - - [13, masure_regression, {t: 0.07, eta: 0.08, kappa: 0.71, delta: -0.8, lambda: 0.41, phi: 0.13}] - - [14, masure_regression, {t: 0.07, eta: 0.07, kappa: 0.69, delta: -0.7, lambda: 0.39, phi: 0.12}] - - [15, masure_regression, {t: 0.07, eta: 0.06, kappa: 0.67, delta: -0.6, lambda: 0.37, phi: 0.11}] - - [16, masure_regression, {t: 0.07, eta: 0.05, kappa: 0.65, delta: -0.5, lambda: 0.35, phi: 0.10}] - - [17, masure_regression, {t: 0.07, eta: 0.04, kappa: 0.63, delta: -0.4, lambda: 0.33, phi: 0.09}] - - [18, masure_regression, {t: 0.07, eta: 0.03, kappa: 0.61, delta: -0.3, lambda: 0.31, phi: 0.08}] - - [19, masure_regression, {t: 0.07, eta: 0.02, kappa: 0.59, delta: -0.2, lambda: 0.29, phi: 0.07}] + - [1, masure_regression, {t: 0.0782, eta: 0.1753, kappa: 0.0955, delta: 0.0, lambda: 0.16, phi: 0.65}] + - [2, masure_regression, {t: 0.0785, eta: 0.1684, kappa: 0.0952, delta: 0.0, lambda: 0.175, phi: 0.65}] + - [3, masure_regression, {t: 0.0769, eta: 0.1776, kappa: 0.0947, delta: 0.0, lambda: 0.145, phi: 0.65}] + - [4, masure_regression, {t: 0.0786, eta: 0.1764, kappa: 0.0933, delta: 0.0, lambda: 0.15, phi: 0.65}] + - [5, masure_regression, {t: 0.0789, eta: 0.1618, kappa: 0.0924, delta: 0.0, lambda: 0.17, phi: 0.65}] + - [6, masure_regression, {t: 0.0788, eta: 0.1756, kappa: 0.0920, delta: 0.0, lambda: 0.155, phi: 0.65}] + - [7, masure_regression, {t: 0.0780, eta: 0.1865, kappa: 0.0894, delta: 0.0, lambda: 0.125, phi: 0.65}] + - [8, masure_regression, {t: 0.0798, eta: 0.1682, kappa: 0.0847, delta: 0.0, lambda: 0.23, phi: 0.65}] + - [9, masure_regression, {t: 0.0804, eta: 0.1801, kappa: 0.0854, delta: 0.0, lambda: 0.15, phi: 0.65}] + - [10, masure_regression, {t: 0.0804, eta: 0.1701, kappa: 0.0834, delta: 0.0, lambda: 0.14, phi: 0.65}] + - [11, masure_regression, {t: 0.0789, eta: 0.1743, kappa: 0.0760, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [12, masure_regression, {t: 0.0842, eta: 0.1976, kappa: 0.0732, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [13, masure_regression, {t: 0.0862, eta: 0.1578, kappa: 0.0667, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [14, masure_regression, {t: 0.0873, eta: 0.2206, kappa: 0.0666, delta: 0.0, lambda: 0.0, phi: 0.65}] + - [15, masure_regression, {t: 0.0883, eta: 0.2052, kappa: 0.0574, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [16, masure_regression, {t: 0.0767, eta: 0.2068, kappa: 0.0549, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [17, masure_regression, {t: 0.0661, eta: 0.1911, kappa: 0.0567, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [18, masure_regression, {t: 0.0614, eta: 0.2122, kappa: 0.0622, delta: 0.0, lambda: 0.10, phi: 0.65}] + - [19, masure_regression, {t: 0.0424, eta: 0.2575, kappa: 0.0677, delta: 0.0, lambda: 0.175, phi: 0.65}] - -bridle_nodes: - # --------------------------------------------------------------- - # headers: - # - id: integer, unique identifier for the node - # - x: x-coordinate [m] - # - y: y-coordinate [m] - # - z: z-coordinate [m] - # - type: node type, either 'knot' or 'pulley' - # --------------------------------------------------------------- - headers: [id, x, y, z, type] - data: - - [21, -0.67653494, 2.10580567, 7.08413599, knot] - - [22, -0.79196791, 1.01882151, 7.57069579, knot] - - [23, -0.79196791, -1.01882151, 7.57069579, knot] - - [24, -0.67653494, -2.10580567, 7.08413599, knot] - - [25, -0.46444477, 0.77505128, 4.23316161, knot] - - [26, -0.46444477, -0.77505128, 4.23316161, knot] - - [27, 0.54344694, 1.57655435, 5.91086207, knot] - - [28, 0.62960591, 0.92728966, 6.64639402, knot] - - [29, 0.62960591, -0.92728966, 6.64639402, knot] - - [30, 0.54344694, -1.57655435, 5.91086207, knot] - - [31, 0.39092886, 0.71221536, 4.26118728, knot] - - [32, 0.39092886, -0.71221536, 4.26118728, knot] - - [33, 0.30679124, 0.49235914, 3.39885672, pulley] - - [34, 0.28625931, 0.0, 3.19030736, knot] - - [35, 0.30679124, -0.49235914, 3.39885672, pulley] - - [36, 0.13269336, 0.23394851, 1.48060245, knot] - - [37, 0.13269336, -0.23394851, 1.48060245, knot] - -bridle_lines: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - rest_length: measured rest length [m] - # - diameter: line diameter [m] - # - material: string, material type (e.g., dyneema) - # - density: material density [kg/m^3] - # --------------------------------------------------------------- - headers: [name, rest_length, diameter, material, density] - data: - # front lines - - [a6, 11.790, 0.002, dyneema,970] - - [A5, 2.870, 0.002, dyneema,970] - - [A3, 3.450, 0.002, dyneema,970] - - [A2, 3.500, 0.002, dyneema,970] - - [A1, 3.690, 0.002, dyneema,970] - - [AII, 3.250, 0.002, dyneema,970] - - [AI, 3.620, 0.002, dyneema,970] - - [amain, 3.910, 0.002, dyneema,970] - # rear lines - - [br1, 4.410, 0.002, dyneema,970] - - [br2, 4.170, 0.002, dyneema,970] - - [br3, 4.140, 0.002, dyneema,970] - - [br5, 3.580, 0.002, dyneema,970] - - [br6, 13.560, 0.002, dyneema,970] - - [BRI, 2.360, 0.002, dyneema,970] - - [BRII, 2.080, 0.002, dyneema,970] - - [BR-main-1, .700, 0.002, dyneema,970] - - [M-line, 2.410, 0.002, dyneema,970] - # - [S, 1.730, 0.002, dyneema,970] # measured - # - [Power Tape, 3.330, 0.002, polyester] # measured - - [Steering Tape, 1.530, 0.002, dyneema,970] # trying 0.2m smalller - - [Power Tape, 2.900, 0.002, dyneema,970] # trying 0.2m smaller - -bridle_connections: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - ci: integer, node id (start) - # - cj: integer, node id (end) - # - ck: integer, third node id (only for pulleys, else omitted or 0) - # --------------------------------------------------------------- - headers: [name, ci, cj, ck] - data: - # front lines - - [a6,1,25,33] #e.g. this is a pulley - - [A5,3,21] - - [A3,5,21] - - [A2,7,22] - - [A1,9,22] - - [A1,11,23] - - [A2,13,23] - - [A3,15,24] - - [A5,17,24] - - [a6,19,26,35] - - [AII,21,25] - - [AI,22,25] - - [AI,23,26] - - [AII,24,26] - - [amain,25,0] - - [amain,26,0] - # rear lines - - [a6,1,33,25] - - [br6,2,33,36] - - [br6,2,36,33] - - [br5,4,27] - - [br3,6,27] - - [br2,8,28] - - [br1,10,28] - - [br1,12,29] - - [br2,14,29] - - [br3,16,30] - - [br5,18,30] - - [a6,19,35,26] - - [br6,20,35,37] - - [br6,20,37,35] - - [BRII,27,31] - - [BRI,28,31] - - [BRI,29,32] - - [BRII,30,32] - - [BR-main-1,31,33] - - [BR-main-1,32,35] - - [M-line,33,36,34] - - [M-line,33,34,36] - - [M-line,35,34,37] - - [M-line,35,37,34] - - [Steering Tape,36,0] - - [Power Tape,34,0] - - [Steering Tape,37,0] diff --git a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_neuralfoil.yaml b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_neuralfoil.yaml similarity index 71% rename from data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_neuralfoil.yaml rename to data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_neuralfoil.yaml index df28c07..171602f 100644 --- a/data/TUDELFT_V3_KITE/CAD_derived_geometry/config_kite_CAD_neuralfoil.yaml +++ b/data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_neuralfoil.yaml @@ -104,126 +104,6 @@ wing_airfoils: -bridle_nodes: - # --------------------------------------------------------------- - # headers: - # - id: integer, unique identifier for the node - # - x: x-coordinate [m] - # - y: y-coordinate [m] - # - z: z-coordinate [m] - # - type: node type, either 'knot' or 'pulley' - # --------------------------------------------------------------- - headers: [id, x, y, z, type] - data: - - [21, -0.67653494, 2.10580567, 7.08413599, knot] - - [22, -0.79196791, 1.01882151, 7.57069579, knot] - - [23, -0.79196791, -1.01882151, 7.57069579, knot] - - [24, -0.67653494, -2.10580567, 7.08413599, knot] - - [25, -0.46444477, 0.77505128, 4.23316161, knot] - - [26, -0.46444477, -0.77505128, 4.23316161, knot] - - [27, 0.54344694, 1.57655435, 5.91086207, knot] - - [28, 0.62960591, 0.92728966, 6.64639402, knot] - - [29, 0.62960591, -0.92728966, 6.64639402, knot] - - [30, 0.54344694, -1.57655435, 5.91086207, knot] - - [31, 0.39092886, 0.71221536, 4.26118728, knot] - - [32, 0.39092886, -0.71221536, 4.26118728, knot] - - [33, 0.30679124, 0.49235914, 3.39885672, pulley] - - [34, 0.28625931, 0.0, 3.19030736, knot] - - [35, 0.30679124, -0.49235914, 3.39885672, pulley] - - [36, 0.13269336, 0.23394851, 1.48060245, knot] - - [37, 0.13269336, -0.23394851, 1.48060245, knot] - -bridle_lines: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - rest_length: measured rest length [m] - # - diameter: line diameter [m] - # - material: string, material type (e.g., dyneema) - # - density: material density [kg/m^3] - # --------------------------------------------------------------- - headers: [name, rest_length, diameter, material, density] - data: - # front lines - - [a6, 11.790, 0.002, dyneema,970] - - [A5, 2.870, 0.002, dyneema,970] - - [A3, 3.450, 0.002, dyneema,970] - - [A2, 3.500, 0.002, dyneema,970] - - [A1, 3.690, 0.002, dyneema,970] - - [AII, 3.250, 0.002, dyneema,970] - - [AI, 3.620, 0.002, dyneema,970] - - [amain, 3.910, 0.002, dyneema,970] - # rear lines - - [br1, 4.410, 0.002, dyneema,970] - - [br2, 4.170, 0.002, dyneema,970] - - [br3, 4.140, 0.002, dyneema,970] - - [br5, 3.580, 0.002, dyneema,970] - - [br6, 13.560, 0.002, dyneema,970] - - [BRI, 2.360, 0.002, dyneema,970] - - [BRII, 2.080, 0.002, dyneema,970] - - [BR-main-1, .700, 0.002, dyneema,970] - - [M-line, 2.410, 0.002, dyneema,970] - # - [S, 1.730, 0.002, dyneema,970] # measured - # - [Power Tape, 3.330, 0.002, polyester] # measured - - [Steering Tape, 1.530, 0.002, dyneema,970] # trying 0.2m smalller - - [Power Tape, 2.900, 0.002, dyneema,970] # trying 0.2m smaller - -bridle_connections: - # --------------------------------------------------------------- - # headers: - # - name: string, line name - # - ci: integer, node id (start) - # - cj: integer, node id (end) - # - ck: integer, third node id (only for pulleys, else omitted or 0) - # --------------------------------------------------------------- - headers: [name, ci, cj, ck] - data: - # front lines - - [a6,1,25,33] #e.g. this is a pulley - - [A5,3,21] - - [A3,5,21] - - [A2,7,22] - - [A1,9,22] - - [A1,11,23] - - [A2,13,23] - - [A3,15,24] - - [A5,17,24] - - [a6,19,26,35] - - [AII,21,25] - - [AI,22,25] - - [AI,23,26] - - [AII,24,26] - - [amain,25,0] - - [amain,26,0] - # rear lines - - [a6,1,33,25] - - [br6,2,33,36] - - [br6,2,36,33] - - [br5,4,27] - - [br3,6,27] - - [br2,8,28] - - [br1,10,28] - - [br1,12,29] - - [br2,14,29] - - [br3,16,30] - - [br5,18,30] - - [a6,19,35,26] - - [br6,20,35,37] - - [br6,20,37,35] - - [BRII,27,31] - - [BRI,28,31] - - [BRI,29,32] - - [BRII,30,32] - - [BR-main-1,31,33] - - [BR-main-1,32,35] - - [M-line,33,36,34] - - [M-line,33,34,36] - - [M-line,35,34,37] - - [M-line,35,37,34] - - [Steering Tape,36,0] - - [Power Tape,34,0] - - [Steering Tape,37,0] - diff --git a/data/TUDELFT_V3_KITE/Deformed_geometry/config_kite_kitesim_deformed.yaml b/data/TUDELFT_V3_KITE/Deformed_geometry/aero_geometry_kitesim_deformed.yaml similarity index 100% rename from data/TUDELFT_V3_KITE/Deformed_geometry/config_kite_kitesim_deformed.yaml rename to data/TUDELFT_V3_KITE/Deformed_geometry/aero_geometry_kitesim_deformed.yaml diff --git a/data/TUDELFT_V3_KITE/Surfplan_derived_geometry/config_kite_surfplan_breukels_regression.yaml b/data/TUDELFT_V3_KITE/Surfplan_derived_geometry/aero_geometry_surfplan_breukels_regression.yaml similarity index 100% rename from data/TUDELFT_V3_KITE/Surfplan_derived_geometry/config_kite_surfplan_breukels_regression.yaml rename to data/TUDELFT_V3_KITE/Surfplan_derived_geometry/aero_geometry_surfplan_breukels_regression.yaml diff --git a/data/TUDELFT_V3_KITE/Surfplan_derived_geometry/config_kite_surfplan_inviscid.yaml b/data/TUDELFT_V3_KITE/Surfplan_derived_geometry/aero_geometry_surfplan_inviscid.yaml similarity index 100% rename from data/TUDELFT_V3_KITE/Surfplan_derived_geometry/config_kite_surfplan_inviscid.yaml rename to data/TUDELFT_V3_KITE/Surfplan_derived_geometry/aero_geometry_surfplan_inviscid.yaml diff --git a/docs/TUDELFT_V3_KITE_reference_frame.png b/docs/TUDELFT_V3_KITE_reference_frame.png new file mode 100644 index 0000000..631a65c Binary files /dev/null and b/docs/TUDELFT_V3_KITE_reference_frame.png differ diff --git a/examples/TUDELFT_V3_KITE/benchmark.py b/examples/TUDELFT_V3_KITE/benchmark.py index ec66ca7..3ba6374 100644 --- a/examples/TUDELFT_V3_KITE/benchmark.py +++ b/examples/TUDELFT_V3_KITE/benchmark.py @@ -83,14 +83,14 @@ def main() -> None: "Breukels regression", data_dir / "CAD_derived_geometry" - / "config_kite_CAD_breukels_regression.yaml", + / "aero_geometry_CAD_breukels_regression.yaml", None, ), ( "Masure ML regression", data_dir / "CAD_derived_geometry" - / "config_kite_CAD_masure_regression.yaml", + / "aero_geometry_CAD_masure_regression.yaml", ml_models_dir, ), ] diff --git a/examples/TUDELFT_V3_KITE/convergence.py b/examples/TUDELFT_V3_KITE/convergence.py index 725ee2e..d8e7baa 100644 --- a/examples/TUDELFT_V3_KITE/convergence.py +++ b/examples/TUDELFT_V3_KITE/convergence.py @@ -16,7 +16,7 @@ def main(): Path(PROJECT_DIR) / "data" / "TUDELFT_V3_KITE" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) convergence_analysis_dir = ( Path(PROJECT_DIR) / "results" / "TUDELFT_V3_KITE" / "convergence_analysis" @@ -84,7 +84,7 @@ def convergence_relation( side_slip=0, yaw_rate=0, spanwise_panel_distribution="uniform", - config_file="config_kite_CAD_CFD_polars.yaml", + config_file="aero_geometry_CAD_CFD_polars.yaml", ): """ Create a 2x3 subplot showing actual values and errors vs n_panels for multiple angles of attack. @@ -448,5 +448,5 @@ def convergence_relation( side_slip=0, yaw_rate=0, spanwise_panel_distribution="uniform", - config_file=Path("CAD_derived_geometry") / "config_kite_CAD_CFD_polars.yaml", + config_file=Path("CAD_derived_geometry") / "aero_geometry_CAD_CFD_polars.yaml", ) diff --git a/examples/TUDELFT_V3_KITE/effect_of_varying_airfoil_shape.py b/examples/TUDELFT_V3_KITE/effect_of_varying_airfoil_shape.py index 7048ea4..0a4597a 100644 --- a/examples/TUDELFT_V3_KITE/effect_of_varying_airfoil_shape.py +++ b/examples/TUDELFT_V3_KITE/effect_of_varying_airfoil_shape.py @@ -37,7 +37,7 @@ def create_parameter_variations(percent_variation=10): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) output_base_dir = ( project_dir / "results" / "TUDELFT_V3_KITE" / "effect_of_varying_airfoil_shape" @@ -161,7 +161,7 @@ def verify_parameter_variations(percent_variation=10): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) variations_base_dir = ( project_dir / "results" / "TUDELFT_V3_KITE" / "effect_of_varying_airfoil_shape" @@ -255,7 +255,7 @@ def generate_2D_plots( project_dir / "data" / "TUDELFT_V3_KITE" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) variations_base_dir = ( project_dir / "results" / "TUDELFT_V3_KITE" / "effect_of_varying_airfoil_shape" @@ -482,7 +482,7 @@ def generate_2D_and_3D_plots( / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) variations_base_dir = ( project_dir / "results" / "TUDELFT_V3_KITE" / "effect_of_varying_airfoil_shape" diff --git a/examples/TUDELFT_V3_KITE/evaluate_stability_derivatives.py b/examples/TUDELFT_V3_KITE/evaluate_stability_derivatives.py index 2a6fec7..332e3ba 100644 --- a/examples/TUDELFT_V3_KITE/evaluate_stability_derivatives.py +++ b/examples/TUDELFT_V3_KITE/evaluate_stability_derivatives.py @@ -3,14 +3,16 @@ from VSM import trim_angle from VSM.core.BodyAerodynamics import BodyAerodynamics from VSM.core.Solver import Solver -from VSM.stability_derivatives import compute_rigid_body_stability_derivatives +from VSM.stability_derivatives import ( + compute_rigid_body_stability_derivatives, + map_derivatives_to_aircraft_frame, +) from VSM.trim_angle import compute_trim_angle -# Default step sizes for finite differences - +# Default step sizes for finite differences PROJECT_DIR = Path(__file__).resolve().parents[2] -n_panels = 30 +n_panels = 100 spanwise_panel_distribution = "uniform" solver_base_version = Solver(reference_point=np.array([0.0, 0.0, 0.0])) @@ -19,7 +21,7 @@ ) body_aero_CAD_CFD_polars = 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_CAD_CFD_polars.yaml"), spanwise_panel_distribution=spanwise_panel_distribution, ) @@ -38,7 +40,6 @@ } # Get reference point from solver for physically correct moment calculations -# The reference point is critical for rotational velocity calculations: # v_rot(r) = omega × (r - r_ref) reference_point = solver_base_version.reference_point @@ -61,7 +62,6 @@ ) print(f"results: {results}") - trim_angle = results["trim_angle"] dCMy_dalpha = results["dCMy_dalpha"] print(f"Trim angle found at {trim_angle:.3f} degrees.") @@ -82,7 +82,7 @@ nondimensionalize_rates=True, # Convert rate derivatives to per hat-rate ) -print("\nComputed stability derivatives:") +print("\nComputed stability derivatives (VSM frame, x rearward, y right, z up):") print("=" * 60) print("Angle derivatives (per radian):") for key, value in derivatives.items(): @@ -93,3 +93,41 @@ for key, value in derivatives.items(): if any(rate in key for rate in ["_dp", "_dq", "_dr"]): print(f" {key}: {value:+.6f}") + +# Apply reference frame transformation using the VSM module function + +derivatives_aircraft = map_derivatives_to_aircraft_frame(derivatives) + +coeffs = ["Cx", "Cy", "Cz", "CMx", "CMy", "CMz"] +angle_keys = [f"d{coeff}_dalpha" for coeff in coeffs] + [ + f"d{coeff}_dbeta" for coeff in coeffs +] +rate_keys = [f"d{coeff}_dp" for coeff in coeffs] + [ + f"d{coeff}_dq" for coeff in coeffs +] + [f"d{coeff}_dr" for coeff in coeffs] + + +def print_combined(title, keys) -> None: + print(f"\n{title}") + header = f" {'Derivative':<18}{'VSM':>14}{'Aircraft':>14}" + print(header) + print(" " + "-" * (len(header) - 2)) + for key in keys: + if key not in derivatives and key not in derivatives_aircraft: + continue + vsm_val = ( + f"{derivatives[key]:+.6f}" if key in derivatives else " n/a" + ) + ac_val = ( + f"{derivatives_aircraft[key]:+.6f}" + if key in derivatives_aircraft + else " n/a" + ) + print(f" {key:<18}{vsm_val:>14}{ac_val:>14}") + + +print_combined("Angle derivatives (per radian)", angle_keys) +print_combined( + "Rate derivatives (per hat-rate, dimensionless)", + rate_keys, +) diff --git a/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py b/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py index 5bbf403..fce951d 100644 --- a/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py +++ b/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py @@ -33,7 +33,7 @@ def main(): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) # bridle_path = ( diff --git a/examples/TUDELFT_V3_KITE/kite_stability_dynamics.py b/examples/TUDELFT_V3_KITE/kite_stability_dynamics.py index e80dff8..6a8bdc0 100644 --- a/examples/TUDELFT_V3_KITE/kite_stability_dynamics.py +++ b/examples/TUDELFT_V3_KITE/kite_stability_dynamics.py @@ -40,7 +40,7 @@ def main(): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) print("=" * 70) diff --git a/examples/TUDELFT_V3_KITE/old/effect_of_tow_point_location.py b/examples/TUDELFT_V3_KITE/old/effect_of_tow_point_location.py deleted file mode 100644 index d8186ff..0000000 --- a/examples/TUDELFT_V3_KITE/old/effect_of_tow_point_location.py +++ /dev/null @@ -1,176 +0,0 @@ -import numpy as np -import time -import matplotlib.pyplot as plt -import pandas as pd -from pathlib import Path -from VSM.core.BodyAerodynamics import BodyAerodynamics -from VSM.core.Solver import Solver -from VSM.fitting import fit_and_evaluate_model -from scipy.interpolate import interp1d - - -def main(): - """ - This script demonstrates how to use the VSM library to perform a 3D aerodynamic analysis of the TUDELFT_V3_KITE. - - The example covers the following steps: - 1. Define file paths for the kite geometry, 2D polars, and bridle geometry. - 2. Load the kite geometry from a CSV file. - 3. Create three BodyAerodynamics objects: - - One using the baseline Breukels input. - - One with corrected polar data. - - One with corrected polar data and bridles. - 4. Initialize the aerodynamic model with a specific wind speed, angle of attack, side slip angle, and yaw rate. - 5. Plot the kite geometry using Matplotlib. - 6. Generate an interactive plot using Plotly. - 7. Plot and save polar curves (both angle of attack and side slip sweeps) for different settings, comparing them to literature data. - """ - - ### 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" - ) - bridle_path = ( - Path(PROJECT_DIR) - / "data" - / "TUDELFT_V3_KITE" - / "CAD_derived_geometry" - / "struc_geometry.yaml" - ) - - ### 2. defining settings - n_panels = 40 - spanwise_panel_distribution = "uniform" - 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, - spanwise_panel_distribution=spanwise_panel_distribution, - bridle_path=bridle_path, - ) - - ### 4. Setting va - Umag = 20 - body_aero_polar_with_bridles.va_initialize(Umag, 8, 5) - - x_tow_point = np.arange(-1.25, 0.5, 0.25) - alpha_range = np.linspace(0, 20, 11) # Angle of attack sweep - ss_range = np.linspace(0, 10, 11) # Sideslip sweep - # Sweep over angle of attack for each reference point and collect all curves - cmz_alpha_all = [] - cmx_alpha_all = [] - cmy_alpha_all = [] - trim_aoa_all = [] - for x in x_tow_point: - solver = Solver(reference_point=[x, 0, 0]) - cmz_alpha = [] - cmx_alpha = [] - cmy_alpha = [] - for alpha in alpha_range: - body_aero_polar_with_bridles.va_initialize(Umag, alpha, 0) - results = solver.solve(body_aero_polar_with_bridles) - cmz_alpha.append(results.get("cmz", np.nan)) - cmx_alpha.append(results.get("cmx", np.nan)) - cmy_alpha.append(results.get("cmy", np.nan)) - # Interpolate to find trim aoa where cmy crosses zero - - cmy_alpha_np = np.array(cmy_alpha) - # Only interpolate if there is a sign change - if np.any(np.diff(np.sign(cmy_alpha_np))): - f_trim = interp1d( - cmy_alpha_np, - alpha_range, - kind="linear", - bounds_error=False, - fill_value="extrapolate", - ) - trim_aoa = f_trim(0.0) - else: - trim_aoa = np.nan - trim_aoa_all.append(trim_aoa) - - cmz_alpha_all.append(cmz_alpha) - cmx_alpha_all.append(cmx_alpha) - cmy_alpha_all.append(cmy_alpha) - - # Triple plot for angle of attack - fig, axs = plt.subplots(1, 3, figsize=(18, 5), sharex=True) - for i, x in enumerate(x_tow_point): - axs[0].plot(alpha_range, cmz_alpha_all[i], label=f"x_tow_point={x:.2f}") - axs[1].plot(alpha_range, cmx_alpha_all[i], label=f"x_tow_point={x:.2f}") - axs[2].plot(alpha_range, cmy_alpha_all[i], label=f"x_tow_point={x:.2f}") - axs[0].set_xlabel("Angle of Attack (deg)") - axs[0].set_ylabel("Yaw Moment Coefficient (cmz)") - axs[0].set_title("cmz vs Angle of Attack") - axs[0].grid() - axs[0].legend() - axs[1].set_xlabel("Angle of Attack (deg)") - axs[1].set_ylabel("Roll Moment Coefficient (cmx)") - axs[1].set_title("cmx vs Angle of Attack") - axs[1].grid() - axs[1].legend() - axs[2].set_xlabel("Angle of Attack (deg)") - axs[2].set_ylabel("Pitch Moment Coefficient (cmy)") - axs[2].set_title("cmy vs Angle of Attack") - axs[2].grid() - axs[2].legend() - plt.suptitle("Dependency on Angle of Attack for Different Reference Points") - plt.tight_layout() - # plt.show() - - # Sweep over sideslip for each reference point and collect all curves - cmz_ss_all = [] - cmx_ss_all = [] - cmy_ss_all = [] - for i, x in enumerate(x_tow_point): - solver = Solver(reference_point=[x, 0, 0]) - cmz_ss = [] - cmx_ss = [] - cmy_ss = [] - for ss in ss_range: - body_aero_polar_with_bridles.va_initialize(Umag, trim_aoa_all[i], ss) - results = solver.solve(body_aero_polar_with_bridles) - cmz_ss.append(results.get("cmz", np.nan)) - cmx_ss.append(results.get("cmx", np.nan)) - cmy_ss.append(results.get("cmy", np.nan)) - cmz_ss_all.append(cmz_ss) - cmx_ss_all.append(cmx_ss) - cmy_ss_all.append(cmy_ss) - - # Triple plot for sideslip - fig, axs = plt.subplots(1, 3, figsize=(18, 5), sharex=True) - for i, x in enumerate(x_tow_point): - axs[0].plot(ss_range, cmz_ss_all[i], label=f"x_tow_point={x:.2f}") - axs[1].plot(ss_range, cmx_ss_all[i], label=f"x_tow_point={x:.2f}") - axs[2].plot(ss_range, cmy_ss_all[i], label=f"x_tow_point={x:.2f}") - axs[0].set_xlabel("Sideslip Angle (deg)") - axs[0].set_ylabel("Yaw Moment Coefficient (cmz)") - axs[0].set_title("cmz vs Sideslip Angle") - axs[0].grid() - axs[0].legend() - axs[1].set_xlabel("Sideslip Angle (deg)") - axs[1].set_ylabel("Roll Moment Coefficient (cmx)") - axs[1].set_title("cmx vs Sideslip Angle") - axs[1].grid() - axs[1].legend() - axs[2].set_xlabel("Sideslip Angle (deg)") - axs[2].set_ylabel("Pitch Moment Coefficient (cmy)") - axs[2].set_title("cmy vs Sideslip Angle") - axs[2].grid() - axs[2].legend() - plt.suptitle("Dependency on Sideslip for Different Reference Points") - plt.tight_layout() - plt.show() - - -if __name__ == "__main__": - main() diff --git a/examples/TUDELFT_V3_KITE/old/tow_angle_vs_aoa.py b/examples/TUDELFT_V3_KITE/old/tow_angle_vs_aoa.py deleted file mode 100644 index dc9d680..0000000 --- a/examples/TUDELFT_V3_KITE/old/tow_angle_vs_aoa.py +++ /dev/null @@ -1,379 +0,0 @@ -import numpy as np -import time -import matplotlib.pyplot as plt -import pandas as pd -from pathlib import Path -from VSM.core.BodyAerodynamics import BodyAerodynamics -from VSM.core.Solver import Solver -from VSM.fitting import fit_and_evaluate_model -from scipy.interpolate import interp1d -from VSM.plot_styling import set_plot_style, PALETTE - -import math - - -def main(): - """ - This script demonstrates how to use the VSM library to perform a 3D aerodynamic analysis of the TUDELFT_V3_KITE. - - The example covers the following steps: - 1. Define file paths for the kite geometry, 2D polars, and bridle geometry. - 2. Load the kite geometry from a CSV file. - 3. Create three BodyAerodynamics objects: - - One using the baseline Breukels input. - - One with corrected polar data. - - One with corrected polar data and bridles. - 4. Initialize the aerodynamic model with a specific wind speed, angle of attack, side slip angle, and yaw rate. - 5. Plot the kite geometry using Matplotlib. - 6. Generate an interactive plot using Plotly. - 7. Plot and save polar curves (both angle of attack and side slip sweeps) for different settings, comparing them to literature data. - """ - - ### 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" - ) - bridle_path = ( - Path(PROJECT_DIR) - / "data" - / "TUDELFT_V3_KITE" - / "CAD_derived_geometry" - / "struc_geometry.yaml" - ) - set_plot_style() - ### 2. defining settings - n_panels = 40 - mass_kite = 12 - spanwise_panel_distribution = "uniform" - reference_point = np.array([0.0, 0, 0]) - solver = Solver(reference_point=reference_point) - - ### 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, - spanwise_panel_distribution=spanwise_panel_distribution, - bridle_path=None, - ) - - ### 4. Setting va - Umag = 20 - yaw_rate = 0 - body_aero_polar_with_bridles.va_initialize(Umag, 5, 0, yaw_rate) - - ### 7. Plotting the polar curves for different angles of attack and side slip angles - # and saving in results with literature - save_folder = Path(PROJECT_DIR) / "results" / "TUDELFT_V3_KITE" - - angle_of_attack_range = np.linspace(-5, 10, 61) - gamma = None - center_of_pressure = np.zeros((len(angle_of_attack_range), 3)) - total_force = np.zeros((len(angle_of_attack_range), 3)) - total_moment = np.zeros((len(angle_of_attack_range), 3)) - cl = np.zeros((len(angle_of_attack_range))) - cd = np.zeros((len(angle_of_attack_range))) - aero_roll = np.zeros((len(angle_of_attack_range))) - tow_angle = np.zeros((len(angle_of_attack_range))) - corner_points = np.array( - [panel.corner_points for panel in body_aero_polar_with_bridles.panels] - ) - begin_time = time.time() - N = len(corner_points) - x_corner = corner_points[N // 2, 0, :] - fbridle_length = np.linalg.norm(x_corner - reference_point) - trim_angle = np.zeros((len(angle_of_attack_range))) - print("Front bridle length:", fbridle_length) - elevation_angle = np.radians(30) - course_angle = np.radians(90) - for i, angle_i in enumerate(angle_of_attack_range): - - body_aero_polar_with_bridles.va_initialize(Umag, angle_i, 0, yaw_rate) - - results = solver.solve(body_aero_polar_with_bridles, gamma_distribution=gamma) - center_of_pressure[i, :] = results["center_of_pressure"] - - x_tow_point = x_corner[0] - center_of_pressure[i, 0] - tow_angle[i] = np.arctan(abs(x_tow_point) / fbridle_length) - - # tow_angle = - total_force[i, :] = np.array( - [ - results["Fx"], - results["Fy"], - results["Fz"], - ] - ) - total_moment[i, :] = np.array( - [ - results["Mx"], - results["My"] - + mass_kite - * 9.81 - * np.cos(elevation_angle) - * np.cos(course_angle) - * 10, - results["Mz"], - ] - ) - - cl[i] = np.sqrt(results["cl"] ** 2 + results["cs"] ** 2) - cd[i] = results["cd"] - aero_roll[i] = np.arctan2(results["cs"], results["cl"]) * 180 / np.pi - - end_time = time.time() - print(f"Time taken for calculations: {end_time - begin_time} seconds") - # Interpolate the angle where My crosses 0 - - # Find indices where My changes sign - sign_change_indices = np.where(np.diff(np.sign(total_moment[:, 1])))[0] - trim_angles = [] - for idx in sign_change_indices: - # Interpolate between angle_of_attack_range[idx] and angle_of_attack_range[idx+1] - f_interp = interp1d( - total_moment[idx : idx + 2, 1], - angle_of_attack_range[idx : idx + 2], - kind="linear", - ) - trim_angle_zero = f_interp(0) - trim_angles.append(trim_angle_zero) - print("Trim angle(s) where My crosses zero:", trim_angles) - dependencies = [ - "np.ones(len(alpha))", - "alpha", - "alpha**2", - ] - # Fit lift coeffcients - fit_cl = fit_and_evaluate_model( - cl, - dependencies=dependencies, - alpha=angle_of_attack_range / 180 * np.pi, - ) - print("Fitted coefficients for lift coefficient:") - print(fit_cl["coeffs"]) - dependencies = [ - "np.ones(len(alpha))", - "alpha", - "alpha**2", - ] - fit_cd = fit_and_evaluate_model( - cd, - dependencies=dependencies, - alpha=angle_of_attack_range / 180 * np.pi, - ) - print("Fitted coefficients for drag coefficient:") - print(fit_cd["coeffs"]) - mean_aoa_exp_ro = 6.37 - std_aoa_exp_ro = 1.93 - mean_aoa_exp_ri = 3 # rough estimate - std_aoa_exp_ri = 0.9 - - plt.figure(figsize=(5, 4)) - plt.plot(angle_of_attack_range, tow_angle * 180 / np.pi) - # Fill between mean ± std for experimental outer and inner - plt.axvspan( - mean_aoa_exp_ro - std_aoa_exp_ro, - mean_aoa_exp_ro + std_aoa_exp_ro, - alpha=0.3, - color=PALETTE["Orange"], - label="Reel-out range", - ) - - plt.axvspan( - mean_aoa_exp_ri - std_aoa_exp_ri, - mean_aoa_exp_ri + std_aoa_exp_ri, - color=PALETTE["Sky Blue"], - alpha=0.3, - label="Reel-in range", - ) - - plt.xlabel(r"Angle of attack, $\alpha_\mathrm{w}$ [$^\circ$]") - plt.ylabel(r"Tow angle, $\lambda_{\mathrm{b}}$ [$^\circ$]") - plt.legend(loc="upper right", frameon=True) - plt.tight_layout() - plt.savefig(save_folder / "tow_angle_vs_aoa.pdf") - plt.show() - # #### 5. Plotting the kite geometry using Matplotlib - # fig,ax = creating_geometry_plot( - # body_aero_polar_with_bridles, - # title="Center of Pressure", - # view_elevation=15, - # view_azimuth=-120, - # ) - # # Plot points in the 3d plot - - # plt.show() - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot(111, projection="3d") - - print("number of panels: ", len(body_aero_polar_with_bridles.panels)) - - for i, panel in enumerate(body_aero_polar_with_bridles.panels): - # Get the corner points of the current panel and close the loop by adding the first point again - x_corners = np.append(corner_points[i, :, 0], corner_points[i, 0, 0]) - y_corners = np.append(corner_points[i, :, 1], corner_points[i, 0, 1]) - z_corners = np.append(corner_points[i, :, 2], corner_points[i, 0, 2]) - - # Plot the panel edges - ax.plot( - x_corners, - y_corners, - z_corners, - color="grey", - label="Panel Edges" if i == 0 else "", - linewidth=1, - ) - for i in range(center_of_pressure.shape[0]): - # center of pressure point - x_cp = center_of_pressure[i, :] - ax.scatter( - x_cp[0], - x_cp[1], - x_cp[2], - s=100, # size of the point - label=f"CP {i+1}", - ) - # Add angle of attack as text near the scatter point - ax.text( - x_cp[0], - x_cp[1] + 0.1, - x_cp[2], - f"{angle_of_attack_range[i]:.1f}°", - color="black", - fontsize=8, - ) - - print("Center of pressure: ", center_of_pressure) - max_diff_x = np.max(center_of_pressure[:, 0]) - np.min(center_of_pressure[:, 0]) - print("max_diff_x: ", max_diff_x / 2.6 * 100, "%") - - # --- Compare designs with different reference point heights --- - ref_heights = np.linspace(0, 8, 8) # Example heights, adjust as needed - # ref_x = np.linspace(0, -0.5, 8) # Example x-coordinates, adjust as needed - - trim_angles = [] - fbridle_lengths = [] - T_half = np.zeros(len(ref_heights)) - T = np.zeros(len(ref_heights)) - t_zero = np.zeros(len(ref_heights)) - - plt.figure() - for i, z in enumerate(ref_heights): - reference_point = np.array([0, 0, z]) - solver = Solver(reference_point=reference_point) - body_aero_polar_with_bridles.va_initialize(Umag, 5, 0, yaw_rate) - corner_points = np.array( - [panel.corner_points for panel in body_aero_polar_with_bridles.panels] - ) - N = len(corner_points) - x_corner = corner_points[N // 2, 0, :] - fbridle_length = np.linalg.norm(x_corner - reference_point) - fbridle_lengths.append(fbridle_length) - - my_vs_alpha = [] - for angle_i in angle_of_attack_range: - body_aero_polar_with_bridles.va_initialize(Umag, angle_i, 0, yaw_rate) - results = solver.solve( - body_aero_polar_with_bridles, gamma_distribution=gamma - ) - my_vs_alpha.append(results["My"]) - dMy_dalpha_field = np.gradient( - my_vs_alpha, angle_of_attack_range - ) # local (central) slope at each grid point - dMy_dalpha = np.interp( - trim_angle[i], angle_of_attack_range, dMy_dalpha_field - ) # slope evaluated at α* - - # α ≈ θ locally ⇒ ∂α/∂θ = 1 ⇒ m_theta = dMy/dα - k = -dMy_dalpha # restoring stiffness [N·m/rad] - print(f"z={z:.2f}, bridle={fbridle_length:.2f}, k={k:.2f}") - - if k <= 0: - T_half[:] = T[:] = t_zero[:] = np.nan # unstable around trim - else: - L = x_corner[2] - z - Iy = mass_kite * (L**2) / 12.0 # slender-rod about CG (adjust if needed) - print(f"Iy={Iy:.2f}") - if Iy <= 0: - T_half[i] = T[i] = t_zero[i] = np.nan - continue - wn = np.sqrt(k / Iy) - T_half[i] = np.pi / wn - T[i] = 2.0 * np.pi / wn - t_zero[i] = np.pi / (2.0 * wn) - print(f"T_half={T_half[i]:.2f}, T={T[i]:.2f}, t_zero={t_zero[i]:.2f}") - plt.plot( - angle_of_attack_range, - my_vs_alpha, - label=f"z={z:.2f}, bridle={fbridle_length:.2f}", - ) - - # Find trim angle (where My crosses zero) - my_vs_alpha = np.array(my_vs_alpha) - sign_change_indices = np.where(np.diff(np.sign(my_vs_alpha)))[0] - if len(sign_change_indices) > 0: - idx = sign_change_indices[0] - f_interp = interp1d( - my_vs_alpha[idx : idx + 2], - angle_of_attack_range[idx : idx + 2], - kind="linear", - ) - trim_angle_zero = f_interp(0) - trim_angles.append(trim_angle_zero) - else: - trim_angles.append(np.nan) - - plt.xlabel("Angle of Attack (degrees)") - plt.ylabel("Moment My") - plt.title("Moment My vs Angle of Attack for Different Reference Point Heights") - plt.legend() - plt.grid() - plt.show() - - # Plot trim angle vs bridle length - plt.figure() - plt.plot(fbridle_lengths, trim_angles, "o-") - plt.xlabel("Bridle Line Length") - plt.ylabel("Trim Angle (deg)") - plt.title("Trim Angle vs Bridle Line Length") - plt.grid() - plt.show() - - fig, axs = plt.subplots(1, 2, figsize=(12, 6)) - axs[0].plot(angle_of_attack_range, cl, label="CL") - axs[0].plot(angle_of_attack_range, fit_cl["data_est"], label="Fitted CL") - axs[0].set_xlabel("Angle of Attack (degrees)") - axs[0].set_ylabel("Lift Coefficient (CL)") - - axs[0].legend() - axs[0].grid() - axs[1].plot(angle_of_attack_range, cd, label="CD") - axs[1].plot(angle_of_attack_range, fit_cd["data_est"], label="Fitted CD") - axs[1].set_xlabel("Angle of Attack (degrees)") - axs[1].set_ylabel("Drag Coefficient (CD)") - - axs[1].legend() - axs[1].grid() - plt.show() - - plt.figure() - plt.plot(fbridle_lengths, T_half, "o-", label="Half Period") - plt.plot(fbridle_lengths, T, "o-", label="Full Period") - # plt.yscale("log") - plt.xlabel("Bridle Length") - plt.ylabel("Period") - plt.title("Period vs Bridle Length") - plt.grid() - plt.legend() - plt.show() - - -if __name__ == "__main__": - main() diff --git a/examples/TUDELFT_V3_KITE/sensivity.py b/examples/TUDELFT_V3_KITE/sensivity.py index dae3c56..53da55b 100644 --- a/examples/TUDELFT_V3_KITE/sensivity.py +++ b/examples/TUDELFT_V3_KITE/sensivity.py @@ -14,7 +14,7 @@ def main(): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) sensitivity_results_dir = ( Path(PROJECT_DIR) / "results" / "TUDELFT_V3_KITE" / "sensitivity_analysis" diff --git a/examples/TUDELFT_V3_KITE/tow_angle_geometry.py b/examples/TUDELFT_V3_KITE/tow_angle_geometry.py index 3809d96..ffa3ae6 100644 --- a/examples/TUDELFT_V3_KITE/tow_angle_geometry.py +++ b/examples/TUDELFT_V3_KITE/tow_angle_geometry.py @@ -36,7 +36,7 @@ def main(): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) save_folder = PROJECT_DIR / "results" / "TUDELFT_V3_KITE" save_folder.mkdir(parents=True, exist_ok=True) diff --git a/examples/TUDELFT_V3_KITE/tow_point_location_parametric_study.py b/examples/TUDELFT_V3_KITE/tow_point_location_parametric_study.py index fad18a7..f7cfac9 100644 --- a/examples/TUDELFT_V3_KITE/tow_point_location_parametric_study.py +++ b/examples/TUDELFT_V3_KITE/tow_point_location_parametric_study.py @@ -36,7 +36,7 @@ def main(): / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" - / "config_kite_CAD_CFD_polars.yaml" + / "aero_geometry_CAD_CFD_polars.yaml" ) bridle_path = ( PROJECT_DIR diff --git a/examples/TUDELFT_V3_KITE/tutorial.py b/examples/TUDELFT_V3_KITE/tutorial.py index c0abddd..783d65f 100644 --- a/examples/TUDELFT_V3_KITE/tutorial.py +++ b/examples/TUDELFT_V3_KITE/tutorial.py @@ -1,4 +1,5 @@ from pathlib import Path +import numpy as np from VSM.core.BodyAerodynamics import BodyAerodynamics from VSM.core.Solver import Solver from VSM.plotting import ( @@ -50,7 +51,7 @@ def main(): ### 2. defining settings n_panels = 50 spanwise_panel_distribution = "uniform" - solver_base_version = Solver() + solver_base_version = Solver(reference_point=np.array([0.0, 0.0, 0.0])) # Step 1: Instantiate BodyAerodynamics objects from different YAML configs cad_derived_geometry_dir = ( @@ -58,25 +59,28 @@ def main(): ) body_aero_CAD_CFD_polars = 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_CAD_CFD_polars.yaml"), spanwise_panel_distribution=spanwise_panel_distribution, ) 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 / "aero_geometry_CAD_CFD_polars.yaml"), spanwise_panel_distribution=spanwise_panel_distribution, - bridle_path=(cad_derived_geometry_dir / "struc_geometry.yaml"), + bridle_path=( + cad_derived_geometry_dir / "struc_geometry_manually_adjusted.yaml" + ), ml_models_dir=(Path(PROJECT_DIR) / "data" / "ml_models"), ) body_aero_CAD_neuralfoil = BodyAerodynamics.instantiate( n_panels=n_panels, - file_path=(cad_derived_geometry_dir / "config_kite_CAD_neuralfoil.yaml"), + file_path=(cad_derived_geometry_dir / "aero_geometry_CAD_neuralfoil.yaml"), spanwise_panel_distribution=spanwise_panel_distribution, ) body_aero_masure_regression = BodyAerodynamics.instantiate( n_panels=n_panels, - file_path=(cad_derived_geometry_dir / "config_kite_CAD_masure_regression.yaml"), + 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, ) diff --git a/examples/machine_learning_for_lei_airfoils/plot_airfoils_in_yaml.py b/examples/machine_learning_for_lei_airfoils/plot_airfoils_in_yaml.py new file mode 100644 index 0000000..cb83d18 --- /dev/null +++ b/examples/machine_learning_for_lei_airfoils/plot_airfoils_in_yaml.py @@ -0,0 +1,319 @@ +#!/usr/bin/env python3 +""" +Plot all airfoils defined in a YAML configuration file. + +This script reads a kite configuration YAML file and plots all the airfoils +defined in the wing_airfoils section. Each airfoil is plotted in a separate +subplot arranged in a single column. +""" + +import yaml +import numpy as np +import matplotlib.pyplot as plt +from pathlib import Path + +from VSM.plot_styling import set_plot_style +from utils_masure_regression_lei_parametric import generate_profile + + +def load_cad_airfoil(dat_file_path): + """ + Load airfoil coordinates from a .dat file. + + Parameters: + ----------- + dat_file_path : Path or str + Path to the .dat file containing airfoil coordinates + + Returns: + -------- + numpy.ndarray or None + Nx2 array of (x, y) coordinates, or None if file doesn't exist + """ + dat_path = Path(dat_file_path) + if not dat_path.exists(): + return None + + try: + # Load CSV data (comma-separated x,y coordinates) + data = np.loadtxt(dat_path, delimiter=",") + return data + except Exception as e: + print(f" Warning: Could not load {dat_path}: {e}") + return None + + +def extract_airfoils_from_yaml(yaml_path): + """ + Extract all airfoil definitions from a YAML configuration file. + + Parameters: + ----------- + yaml_path : Path or str + Path to the YAML configuration file + + Returns: + -------- + list + List of dictionaries, each containing airfoil parameters + """ + with open(yaml_path, "r", encoding="utf-8") as f: + config = yaml.safe_load(f) + + airfoils = [] + + if "wing_airfoils" in config and "data" in config["wing_airfoils"]: + wing_airfoils_data = config["wing_airfoils"]["data"] + + for airfoil_entry in wing_airfoils_data: + airfoil_id, airfoil_type, params = airfoil_entry + + if airfoil_type == "masure_regression": + airfoil_info = { + "id": airfoil_id, + "type": airfoil_type, + "t": params.get("t", 0.1), + "eta": params.get("eta", 0.3), + "kappa": params.get("kappa", 0.1), + "delta": params.get("delta", 0.0), + "lambda": params.get("lambda", 0.5), + "phi": params.get("phi", 0.7), + } + airfoils.append(airfoil_info) + + return airfoils + + +def plot_all_airfoils(yaml_path, output_path=None, cad_airfoils_dir=None): + """ + Plot all airfoils from a YAML configuration file. + + Parameters: + ----------- + yaml_path : Path or str + Path to the YAML configuration file + output_path : Path or str, optional + Path to save the output figure + cad_airfoils_dir : Path or str, optional + Directory containing CAD-sliced airfoil .dat files (named 1.dat, 2.dat, etc.) + """ + # Extract airfoils from YAML + airfoils = extract_airfoils_from_yaml(yaml_path) + + if not airfoils: + print("No airfoils found in the YAML file.") + return + + n_airfoils = len(airfoils) + print(f"Found {n_airfoils} airfoil(s) in {yaml_path}") + + # Create figure with subplots + fig_width = 12 + fig_height = 2 * n_airfoils + fig, axes = plt.subplots(n_airfoils, 1, figsize=(fig_width, fig_height)) + + # Handle single airfoil case + if n_airfoils == 1: + axes = [axes] + + # Plot each airfoil + for idx, (ax, airfoil) in enumerate(zip(axes, airfoils)): + print(f"\nGenerating airfoil {idx + 1}/{n_airfoils}:") + print(f" ID: {airfoil['id']}") + print( + f" Parameters: t={airfoil['t']:.4f}, η={airfoil['eta']:.4f}, " + f"κ={airfoil['kappa']:.4f}, δ={airfoil['delta']:.1f}°, " + f"λ={airfoil['lambda']:.4f}, φ={airfoil['phi']:.4f}" + ) + + try: + # Generate airfoil geometry using complete implementation + result = generate_profile( + t_val=airfoil["t"], + eta_val=airfoil["eta"], + kappa_val=airfoil["kappa"], + delta_val=airfoil["delta"], + lambda_val=airfoil["lambda"], + phi_val=airfoil["phi"], + ) + + # Unpack result (returns tuple: all_points, profile_name, seam_a) + if result is None or len(result) != 3: + raise ValueError("generate_profile returned invalid result") + + points, _, _ = result + + # Plot generated airfoil (black) + ax.plot( + points[:, 0], + points[:, 1], + "k-", + linewidth=2, + label="Generated (masure_regression)", + ) + ax.fill(points[:, 0], points[:, 1], color="lightgray", alpha=0.3) + + # Try to load and plot CAD airfoil if available (blue) + if cad_airfoils_dir is not None: + cad_file = Path(cad_airfoils_dir) / f"{airfoil['id']}.dat" + cad_points = load_cad_airfoil(cad_file) + + if cad_points is not None: + ax.plot( + cad_points[:, 0], + cad_points[:, 1], + "b-", + linewidth=2, + label="CAD (sliced)", + alpha=0.8, + ) + print(f" ✓ CAD airfoil loaded from {cad_file.name}") + max_y = max(np.concatenate((points[:, 1], cad_points[:, 1]))) + min_y = min(np.concatenate((points[:, 1], cad_points[:, 1]))) + else: + max_y = max(points[:, 1]) + min_y = min(points[:, 1]) + # Formatting + ax.set_aspect("equal", "box") + ax.grid(True, alpha=0.3) + ax.set_xlabel("$x/c$ (-)") + ax.set_ylabel("$y/c$ (-)") + + # Add legend if CAD data was plotted + if cad_airfoils_dir is not None: + ax.legend(loc="upper right", fontsize=8) + + # Title with parameters + title = ( + f"Airfoil {airfoil['id']}: " + f"$t$={airfoil['t']:.3f}, $\\eta$={airfoil['eta']:.3f}, " + f"$\\kappa$={airfoil['kappa']:.3f}, $\\delta$={airfoil['delta']:.1f}°, " + f"$\\lambda$={airfoil['lambda']:.3f}, $\\phi$={airfoil['phi']:.3f}" + ) + ax.set_title(title, pad=10) + + # Set reasonable axis limits + y_range = max_y - min_y + ax.set_xlim(-0.05, 1.05) + ax.set_ylim(min_y - 0.1 * y_range, max_y + 0.1 * y_range) + + print(f" ✓ Successfully generated") + + except Exception as e: + print(f" ✗ Error generating airfoil: {e}") + ax.text( + 0.5, + 0.5, + f"Error: {str(e)}", + ha="center", + va="center", + transform=ax.transAxes, + ) + ax.set_aspect("equal", "box") + ax.grid(True, alpha=0.3) + + # Adjust layout + plt.tight_layout() + + # Save figure if output path provided + if output_path: + output_path = Path(output_path) + output_path.parent.mkdir(parents=True, exist_ok=True) + plt.tight_layout() + plt.savefig(output_path) + print(f"\nFigure saved to: {output_path}") + + +def main(): + """ + Main function to plot airfoils from a YAML configuration file. + + Modify the yaml_path variable below to point to your desired YAML file. + """ + # Apply plot styling + set_plot_style() + PROJECT_DIR = Path(__file__).resolve().parents[2] + + # ===== USER CONFIGURATION ===== + # Modify this path to your YAML file + yaml_path = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "CAD_derived_geometry" + / "aero_geometry_CAD_masure_regression.yaml" + ) + + # Optional: specify output path for saving the figure + output_path = ( + Path(PROJECT_DIR) + / "results" + / "TUDELFT_V3_KITE" + / f"airfoils_in_{yaml_path.stem}.pdf" + ) + + # Optional: directory containing CAD-sliced airfoils (1.dat, 2.dat, etc.) + cad_airfoils_dir = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "CAD_derived_geometry" + / "airfoils_sliced_from_CAD" + ) + # ============================== + + print("=" * 60) + print("Plotting all airfoils from YAML configuration") + print("=" * 60) + print(f"YAML file: {yaml_path}") + + if not yaml_path.exists(): + print(f"ERROR: YAML file not found at {yaml_path}") + return + + # First, let's print what's in the YAML to debug + print("\nReading YAML file...") + airfoils = extract_airfoils_from_yaml(yaml_path) + + if airfoils: + print(f"\nFound {len(airfoils)} airfoil(s):") + for airfoil in airfoils: + print(f"\n Airfoil ID {airfoil['id']}:") + print(f" t (tube_size) = {airfoil['t']:.6f}") + print(f" eta (c_x, camber position) = {airfoil['eta']:.6f}") + print(f" kappa (c_y, camber height) = {airfoil['kappa']:.6f}") + print(f" delta (TE_angle) = {airfoil['delta']:.6f}") + print(f" lambda (TE_cam_tension) = {airfoil['lambda']:.6f}") + print(f" phi (LE_tension) = {airfoil['phi']:.6f}") + + # Check for unreasonable values + if airfoil["kappa"] > 0.3: + print( + f" WARNING: kappa={airfoil['kappa']:.6f} is very large! Should be < 0.3" + ) + if airfoil["eta"] > 0.5: + print( + f" WARNING: eta={airfoil['eta']:.6f} is very large! Should be < 0.5" + ) + if airfoil["t"] > 0.2: + print( + f" WARNING: t={airfoil['t']:.6f} is very large! Should be < 0.2" + ) + + # Check if CAD airfoils directory exists + if cad_airfoils_dir.exists(): + print(f"\nCAD airfoils directory found: {cad_airfoils_dir}") + else: + print(f"\nWarning: CAD airfoils directory not found: {cad_airfoils_dir}") + cad_airfoils_dir = None + + # Generate plots + plot_all_airfoils(yaml_path, output_path, cad_airfoils_dir) + + print("\n" + "=" * 60) + print("Done!") + print("=" * 60) + + +if __name__ == "__main__": + main() diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index 9a45f45..e8051f6 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -139,7 +139,7 @@ def instantiate( Returns: BodyAerodynamics instance. - YAML file structure expected (see config_kite_surfplan.yaml): + YAML file structure expected (see aero_geometry_surfplan_inviscid.yaml): - wing_sections: headers: [airfoil_id, LE_x, LE_y, LE_z, TE_x, TE_y, TE_z] @@ -1329,7 +1329,7 @@ def va_initialize( Parameters: Umag (float): Magnitude of the velocity. angle_of_attack (float): Angle of attack in degrees. - side_slip (float): Sideslip angle in degrees, a minus is added because its defined counter-clockwise. + side_slip (float): Sideslip angle in degrees. yaw_rate (float): Yaw rate about the body z-axis, default is 0.0. pitch_rate (float): Pitch rate about the body y-axis, default is 0.0. roll_rate (float): Roll rate about the body x-axis, default is 0.0. @@ -1337,8 +1337,7 @@ def va_initialize( """ # Convert angles to radians aoa_rad = np.deg2rad(angle_of_attack) - # a (-) is added because its defined counter-clockwise - side_slip_rad = -np.deg2rad(side_slip) + side_slip_rad = np.deg2rad(side_slip) # Set the va attribute using the setter with keyword arguments # The va setter requires keyword-only arguments after the va parameter diff --git a/src/VSM/plotting.py b/src/VSM/plotting.py index 5de85c9..e2ef621 100644 --- a/src/VSM/plotting.py +++ b/src/VSM/plotting.py @@ -399,7 +399,7 @@ def generate_3D_polar_data( begin_time = time.time() results = solver.solve(body_aero, gamma_distribution=gamma) run_time[i] = time.time() - begin_time - print(f"Angle: {angle_i:.1f}deg. Time: {run_time[i]:.4f}s") + print(f"{angle_type}: {angle_i:.1f}deg. Time: {run_time[i]:.4f}s") cl[i] = results["cl"] cd[i] = results["cd"] cs[i] = results["cs"] @@ -629,8 +629,9 @@ def plot_polars( if not os.path.exists(save_path): os.makedirs(save_path) # Save the plot with the specified title and data type - plt.savefig(Path(save_path) / f"{title}{data_type}", bbox_inches="tight") - print(f"Plot saved as: {save_path}") + save_path_ = Path(save_path) / f"{title}{data_type}" + plt.savefig(save_path_, bbox_inches="tight") + print(f"Plot saved as: {save_path_}") # showing plot if is_show: diff --git a/src/VSM/stability_derivatives.py b/src/VSM/stability_derivatives.py index f2fa192..8fbe1e4 100644 --- a/src/VSM/stability_derivatives.py +++ b/src/VSM/stability_derivatives.py @@ -43,7 +43,8 @@ def compute_rigid_body_stability_derivatives( angle_of_attack : float Baseline angle of attack in degrees. side_slip : float - Baseline sideslip angle in degrees (positive for starboard-to-port flow). + 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 @@ -112,8 +113,8 @@ def compute_rigid_body_stability_derivatives( param_names = ("alpha", "beta", "p", "q", "r") default_steps = { - "alpha": 0.5, # degrees - "beta": 0.5, # degrees + "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 @@ -296,3 +297,80 @@ def _central_difference( 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 diff --git a/tests/AirfoilAerodynamics/test_airfoil_aerodynamics.py b/tests/AirfoilAerodynamics/test_airfoil_aerodynamics.py index 0e1a688..38d0a1d 100644 --- a/tests/AirfoilAerodynamics/test_airfoil_aerodynamics.py +++ b/tests/AirfoilAerodynamics/test_airfoil_aerodynamics.py @@ -31,7 +31,7 @@ def masure_test_config(): """Configuration path for masure_regression tests.""" test_dir = Path(__file__).parent.parent.parent return str( - test_dir / "data" / "TUDELFT_V3_KITE" / "config_kite_CAD_masure_regression.yaml" + test_dir / "data" / "TUDELFT_V3_KITE" / "aero_geometry_CAD_masure_regression.yaml" ) diff --git a/tests/AirfoilAerodynamics/test_cache_cleanup.py b/tests/AirfoilAerodynamics/test_cache_cleanup.py index e0b7497..892085e 100644 --- a/tests/AirfoilAerodynamics/test_cache_cleanup.py +++ b/tests/AirfoilAerodynamics/test_cache_cleanup.py @@ -19,7 +19,7 @@ def config_path(self): test_dir / "data" / "TUDELFT_V3_KITE" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) @pytest.fixture diff --git a/tests/AirfoilAerodynamics/test_masure.py b/tests/AirfoilAerodynamics/test_masure.py index fd0eda5..76abfc7 100644 --- a/tests/AirfoilAerodynamics/test_masure.py +++ b/tests/AirfoilAerodynamics/test_masure.py @@ -18,7 +18,7 @@ def masure_config_path(self): test_dir / "data" / "TUDELFT_V3_KITE" - / "config_kite_CAD_masure_regression.yaml" + / "aero_geometry_CAD_masure_regression.yaml" ) @pytest.fixture diff --git a/tests/BodyAerodynamics/config_kite_CAD_inviscid.yaml b/tests/BodyAerodynamics/aero_geometry_CAD_inviscid.yaml similarity index 100% rename from tests/BodyAerodynamics/config_kite_CAD_inviscid.yaml rename to tests/BodyAerodynamics/aero_geometry_CAD_inviscid.yaml diff --git a/tests/BodyAerodynamics/test_kite_global_properties.py b/tests/BodyAerodynamics/test_kite_global_properties.py index 1bf3e99..b98064f 100644 --- a/tests/BodyAerodynamics/test_kite_global_properties.py +++ b/tests/BodyAerodynamics/test_kite_global_properties.py @@ -19,7 +19,7 @@ class TestKiteGlobalProperties: def config_file_path(self): """Path to the test configuration file.""" test_dir = Path(__file__).parent - return test_dir / "config_kite_CAD_inviscid.yaml" + return test_dir / "aero_geometry_CAD_inviscid.yaml" @pytest.fixture def body_aero(self, config_file_path): @@ -190,7 +190,7 @@ class TestGeometricPropertiesWithDifferentPanelCounts: def config_file_path(self): """Path to the test configuration file.""" test_dir = Path(__file__).parent - return test_dir / "config_kite_CAD_inviscid.yaml" + return test_dir / "aero_geometry_CAD_inviscid.yaml" @pytest.mark.parametrize("n_panels", [20, 30, 40, 50]) def test_geometric_properties_mesh_independent(self, config_file_path, n_panels):