diff --git a/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Re1e6_beta_sweep_alpha_13_Vire2022_CorrectedByPoland2025.csv b/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Re1e6_beta_sweep_alpha_13_Vire2022_CorrectedByPoland2025.csv new file mode 100644 index 0000000..29a30ff --- /dev/null +++ b/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Re1e6_beta_sweep_alpha_13_Vire2022_CorrectedByPoland2025.csv @@ -0,0 +1,5 @@ +alpha,beta,CL,CD,CS +13.02,0.0,1.049175234495937,0.1108439499820641,2.5243354512823105e-05 +13.02,4.0,1.0343396181060318,0.1188194181264538,0.040520130466239486 +13.02,8.0,0.9941161938410832,0.1412429431126459,0.07897220142516528 +13.02,12.0,0.7945103258966532,0.2292754151113001,0.07171027635776697 diff --git a/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Rey_10e5_Poland2025_beta_sweep_alpha_13_02.csv b/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Rey_10e5_Poland2025_beta_sweep_alpha_13_02.csv deleted file mode 100644 index b25bc9c..0000000 --- a/data/TUDELFT_V3_KITE/3D_polars_literature/CFD_RANS_Rey_10e5_Poland2025_beta_sweep_alpha_13_02.csv +++ /dev/null @@ -1,5 +0,0 @@ -alpha,beta,CL,CD,CS,CMx,CMy,CMz -13.02,0,1.049472106225422,0.1109605616360219,2.4385142992888863e-05,0.0002462075273015,-0.3256966363266202,-0.000111600945403 -13.02,4,1.034347132862086,0.1187231349047993,-0.040581606940825916,-0.0989628231010319,-0.3252096762447942,-0.0635164012620871 -13.02,8,0.9946801617025854,0.1413864496187794,-0.07900689629021146,-0.1933768769462073,-0.3261046083196199,-0.1265382345550858 -13.02,12,0.9328481343698932,0.1769306090870539,-0.11428499995717821,-0.2874548226069351,-0.3295173645788699,-0.1855631093369855 diff --git a/data/TUDELFT_V3_KITE/3D_polars_literature/V3_CL_CD_CS_alpha_sweep_for_beta_0_WindTunnel_Poland_2025_Rey_560e4.csv b/data/TUDELFT_V3_KITE/3D_polars_literature/V3_CL_CD_CS_alpha_sweep_for_beta_0_WindTunnel_Poland_2025_Rey_560e4.csv deleted file mode 100644 index cfdacb7..0000000 --- a/data/TUDELFT_V3_KITE/3D_polars_literature/V3_CL_CD_CS_alpha_sweep_for_beta_0_WindTunnel_Poland_2025_Rey_560e4.csv +++ /dev/null @@ -1,18 +0,0 @@ -CL,CL_ci,CD,CD_ci,alpha --0.280233153560818,0.0110705335388759,0.0978529777593977,0.0129358291967856,-11.5681211570017 --0.214491373614373,0.0122872270418091,0.0685357795740722,0.0155251384847079,-6.09905950161158 --0.00076803912330013,0.0106850856430888,0.0488586773900786,0.00766964551596273,-1.99963855771642 -0.0739899513915977,0.00953159417113588,0.0573908146562723,0.00888991109005732,-1.33481996708469 -0.465192674677979,0.00947570706047279,0.0561678776611973,0.00908481079346789,3.08107846652584 -0.610180915884597,0.00452502283166885,0.0743102433724348,0.00462168981528609,5.41284642026105 -0.743035635092799,0.00946575579590553,0.0866430637231721,0.00963341345922823,7.35032445798279 -0.887298865863501,0.0106380475782356,0.102336952277841,0.00900521264610847,9.38243360452917 -0.925565515163773,0.0142384313990637,0.175994396308134,0.0150712042894669,11.4644251663019 -0.932725198798349,0.0158775060283244,0.188886748057569,0.0169368686663311,12.4610557905447 -0.951228281361335,0.014573344990979,0.21274544018436,0.0180470531864394,13.3523481658782 -0.977071520361125,0.0181030479414339,0.227886839663643,0.0204637541133525,14.540186234232 -1.00856267581916,0.0166588473137047,0.24447585438451,0.0183782349712973,16.2253663705088 -1.06810346425319,0.0154900038338791,0.319919197381967,0.0229888926456552,18.2973462373086 -1.00999380129617,0.0154476585577098,0.35730163245484,0.0151447117903846,20.2246928771348 -0.997963560484868,0.0118111212863406,0.419985118207109,0.0116052986279285,23.0303543565816 -0.974927517924153,0.0130172960362764,0.434984255224326,0.0121749295096374,24.5411952103548 diff --git a/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_alpha_sweep_beta_0_Poland2025.csv b/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_alpha_sweep_beta_0_Poland2025.csv new file mode 100644 index 0000000..3805894 --- /dev/null +++ b/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_alpha_sweep_beta_0_Poland2025.csv @@ -0,0 +1,18 @@ +alpha,beta,CL,CL_ci,CD,CD_ci,CMx,CMx_ci,CMy,CMy_ci,CMz,CMz_ci +-11.568121157001665,0,-0.28023315356081796,0.011070533538875897,0.09785297775939765,0.012935829196785585,0.00437157475590111,0.0286993060001136,0.0185274596512998,0.0282608431382968,0.0210669698944731,0.0271779286988212 +-6.099059501611582,0,-0.2144913736143728,0.012287227041809052,0.06853577957407217,0.015525138484707877,-0.00755498573046038,0.0360392337568146,0.0114795207949488,0.032500454158528,0.000755708546556479,0.0325068410211215 +-1.9996385577164184,0,-0.0007680391233001299,0.01068508564308875,0.04885867739007862,0.0076696455159627275,-0.0158988564524073,0.037171452190456,-0.00330931153035509,0.0341996655373784,-0.00759575821083012,0.0318738692894563 +-1.3348199670846914,0,0.07398995139159774,0.00953159417113588,0.05739081465627228,0.008889911090057321,-0.0100130430688713,0.0325822709227836,-0.040806162473867,0.0288473702788264,0.0123645037540894,0.0280588394174328 +3.0810784665258444,0,0.46519267467797903,0.009475707060472788,0.05616787766119728,0.009084810793467887,-0.0258618167462268,0.0224844541725969,0.0308738607806185,0.0284283853889457,0.0140251297074127,0.0210561286933184 +5.412846420261046,0,0.610180915884597,0.0045250228316688484,0.07431024337243476,0.004621689815286087,-0.0221407525325612,0.0156450069084061,0.0144224945436694,0.0145159160955344,0.0229736195850115,0.0131474110678698 +7.350324457982788,0,0.7430356350927991,0.009465755795905533,0.08664306372317211,0.009633413459228232,-0.0312824340278635,0.0266562268673516,0.0563725541445889,0.0308395985679497,0.0178991263350178,0.0229903897494152 +9.382433604529174,0,0.887298865863501,0.010638047578235603,0.10233695227784106,0.00900521264610847,-0.02336356591712,0.030224035077201,0.0392859550909577,0.0338451478858252,0.0251841910742331,0.0243551218371345 +11.464425166301869,0,0.9255655151637727,0.014238431399063743,0.17599439630813363,0.01507120428946693,-0.0228070181352731,0.0334609894219685,0.0315587237644146,0.0416478984370718,0.0166303297008157,0.0280255456277692 +12.461055790544702,0,0.9327251987983486,0.015877506028324398,0.18888674805756867,0.01693686866633114,-0.00995600648602306,0.0315207782517628,0.075218719355394,0.0445381386567502,0.00622294840687863,0.0280185877752834 +13.352348165878231,0,0.9512282813613351,0.014573344990978988,0.2127454401843595,0.018047053186439406,-0.0114210223486002,0.0312159263852684,0.0258594620554227,0.0436598327054732,0.0184960316024415,0.0273096754789724 +14.540186234231973,0,0.9770715203611249,0.018103047941433924,0.22788683966364287,0.0204637541133525,-0.0155344068837501,0.032163181237033,0.0531223656310487,0.0520680554723465,0.00865439721248052,0.0284861083051677 +16.225366370508798,0,1.008562675819161,0.016658847313704747,0.24447585438450992,0.018378234971297258,-0.00039395697021092,0.027685485343775,0.0886930289387169,0.0518597446938309,0.00169172965990106,0.0253734777646321 +18.29734623730859,0,1.0681034642531901,0.0154900038338791,0.319919197381967,0.02298889264565519,-0.00860452331661369,0.0325209573901093,0.0403243136368007,0.0469956298727439,0.0390052407209399,0.0289388526822307 +20.224692877134817,0,1.0099938012961718,0.015447658557709807,0.3573016324548397,0.015144711790384626,-0.014479896219406,0.0355487235483263,0.0941951402601613,0.0502408375714964,0.00581594862390604,0.0314948818234675 +23.03035435658158,0,0.9979635604848679,0.01181112128634061,0.4199851182071093,0.011605298627928523,-0.0316878321298153,0.0342147018358101,0.0428177136177952,0.038177499309995,0.000520398672799248,0.0403006191993008 +24.54119521035482,0,0.9749275179241531,0.01301729603627638,0.43498425522432616,0.012174929509637398,-0.0274762633692519,0.0307244717194419,0.0610071490706344,0.0432779678778292,0.00446788426356059,0.0247795303139092 diff --git a/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_beta_sweep_alpha_13_Poland2025.csv b/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_beta_sweep_alpha_13_Poland2025.csv new file mode 100644 index 0000000..a3b8c54 --- /dev/null +++ b/data/TUDELFT_V3_KITE/3D_polars_literature/WindTunnel_Re5e5_beta_sweep_alpha_13_Poland2025.csv @@ -0,0 +1,18 @@ +alpha,beta,CL,CL_ci,CD,CD_ci,CS,CS_ci,CMx,CMx_ci,CMy,CMy_ci,CMz,CMz_ci +12.5,-20.03306181832601,0.777951405622742,0.008796856833471997,0.20965299417515762,0.007808202513084929,-0.07250478256562636,0.008217489683777739,-0.142672196724445,0.0239788146347902,0.105497635459822,0.0267330242074815,0.170095068206459,0.0225675769240559 +12.5,-14.015708749779463,0.8302403074977395,0.012955285139438909,0.1931747523114195,0.009309425805962184,-0.034449390408806985,0.010246626789327874,-0.0715148071806949,0.0245170436517844,0.150692119741533,0.0368042775342016,0.135075597754474,0.0256745216234829 +12.5,-12.01532688930784,0.8486356648182445,0.011645372133727057,0.1845659041133973,0.009763266656747276,-0.033611967911577785,0.00813554565837993,-0.0629870090871866,0.0231105680939916,0.15437286560363,0.0329059321452925,0.10643956020603,0.0234338941148544 +12.5,-10.02193052375522,0.844372693044939,0.011365030678061794,0.18101829425148508,0.0106356859804954,-0.048093781193258786,0.008553952752862554,-0.0634736057356462,0.0242190760151681,0.157390562833406,0.0329547263093356,0.0566738411783359,0.0236325973188328 +12.5,-8.055560834151049,1.006523668797301,0.010837935241264824,0.10802947323957199,0.010272790394690522,-0.12184527056447772,0.008895700961960736,-0.11691206948282,0.0267860169947592,0.180125302822093,0.0348793986829733,0.01411451492199,0.0224436797592357 +12.5,-6.029862930104683,0.9537384196530716,0.015370702522261782,0.1599950502942629,0.021244107617413247,-0.06548959989623337,0.010741698337113218,-0.0408294278667635,0.0286263704955798,0.130546219499605,0.0402005819509749,0.0427415331507237,0.0285571632407517 +12.5,-4.002438813376063,0.9184046713894765,0.014982732979344607,0.17243465965805976,0.01769809717193447,-0.005348333591515309,0.010089797365083922,0.0115795457849802,0.0321603795235286,0.126773298197394,0.0393885813343565,0.0844651336782557,0.0290144489338355 +12.5,-1.9904369199291798,0.9186803652592959,0.014083307882746029,0.18249542351422535,0.015117905473731102,0.020971896776982565,0.011043933305536349,0.0289340426385215,0.0329564536685998,0.111899049073512,0.0421249007013808,0.0701542646332426,0.0286118956952176 +12.5,0.006236934639288943,0.9327251987983486,0.015877506028324398,0.18888674805756867,0.01693686866633114,0.01367763821815797,0.0104799989282744,0.00995600648602306,0.0315207782517628,0.075218719355394,0.0445381386567502,-0.00622294840687863,0.0280185877752834 +12.5,2.0092166497800514,0.9271547044901232,0.015478908445309329,0.18918201422866804,0.018493063677438935,0.020212172896745345,0.010114147444868611,0.00380194605735361,0.030972772935572,0.0670978249130965,0.0424632374852895,-0.0753837021818302,0.0278145126282776 +12.5,4.026042853894021,0.9264779163925219,0.014463729337635955,0.19053897242035212,0.015969671394721947,0.05711214792710673,0.00981487690616588,0.0263085549643294,0.0281338870066655,0.0775354510272678,0.0442854459217499,-0.0721693066715767,0.0275663442296409 +12.5,6.032313484777512,0.904332198504447,0.013675018771545895,0.19875546429064425,0.015903775687221792,0.07086368222790036,0.009670037966489518,0.0462877892104699,0.0284111602498085,0.108084402998637,0.0423601024915454,-0.0755187863584942,0.0249358041317863 +12.5,8.040287349192912,0.866938735025727,0.013157295125787564,0.1832961756463486,0.014969734589144206,0.08835041873904856,0.009641887208934047,0.0889167401541545,0.028995113913509,0.145946798410035,0.0393553788984878,-0.0454755154487997,0.0282073616723269 +12.5,10.037892784949944,0.8939409571719724,0.01295305280419159,0.1904967024240417,0.013890898029644879,0.08309912378412976,0.01027379051446907,0.0940849815137015,0.0260604537199178,0.124385920940314,0.0361359218663154,-0.0706091247819989,0.0277191625554416 +12.5,12.012352613272146,0.8654661435745437,0.013292408593805123,0.1951719602829853,0.014118408150101295,0.027089361225773816,0.010568941682094995,0.0369289869110848,0.0289402501888588,0.142488643956464,0.0388826159056847,-0.165668764638255,0.0296007824043359 +12.5,14.03755227969408,0.8770062636432303,0.01279498307188091,0.18972929247964704,0.013000904667972327,0.0823523935439637,0.010195284111974946,0.103505312128987,0.0250040737685825,0.116472214689101,0.0383725081672792,-0.122119303392468,0.0295713217237066 +12.5,20.050484500008725,0.839104780120975,0.010512431660821126,0.2263726395107008,0.009036578817891099,0.11071283678269149,0.00811056514566816,0.142420967084123,0.0233827963384127,0.0117062004525746,0.0325559862184502,-0.142771430598793,0.022400434567836 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/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() diff --git a/examples/TUDELFT_V3_KITE/benchmark_runtime.py b/examples/TUDELFT_V3_KITE/benchmark_runtime.py new file mode 100644 index 0000000..3ba6374 --- /dev/null +++ b/examples/TUDELFT_V3_KITE/benchmark_runtime.py @@ -0,0 +1,174 @@ +import time +from pathlib import Path +from typing import Dict, List, Tuple + +from VSM.core.BodyAerodynamics import BodyAerodynamics +from VSM.core.Solver import Solver + + +def instantiate_body( + config_path: Path, + n_panels: int, + spanwise_panel_distribution: str, + ml_models_dir: Path | None = None, +) -> BodyAerodynamics: + if not config_path.exists(): + raise FileNotFoundError(f"Configuration file not found: {config_path}") + if ml_models_dir is not None and not ml_models_dir.exists(): + raise FileNotFoundError(f"ml_models_dir does not exist: {ml_models_dir}") + + kwargs: Dict[str, Path] = {} + if ml_models_dir is not None: + kwargs["ml_models_dir"] = ml_models_dir + + return BodyAerodynamics.instantiate( + n_panels=n_panels, + file_path=config_path, + spanwise_panel_distribution=spanwise_panel_distribution, + **kwargs, + ) + + +def run_sweep( + body_aero: BodyAerodynamics, + solver: Solver, + Umag: float, + alpha_values: List[float], + side_slip: float, + yaw_rate: float, + warm_start: bool, +) -> Tuple[List[float], List[Dict[str, float]]]: + gamma_prev = None + timings: List[float] = [] + results: List[Dict[str, float]] = [] + + for alpha in alpha_values: + body_aero.va_initialize(Umag, alpha, side_slip, yaw_rate) + start_time = time.perf_counter() + if warm_start: + solver_results = solver.solve(body_aero, gamma_distribution=gamma_prev) + else: + solver_results = solver.solve(body_aero) + elapsed = time.perf_counter() - start_time + timings.append(elapsed) + results.append( + { + "alpha": alpha, + "cl": solver_results.get("cl"), + "cd": solver_results.get("cd"), + "cs": solver_results.get("cs"), + } + ) + gamma_prev = solver_results.get("gamma_distribution") if warm_start else None + + return timings, results + + +def summarise_timings(label: str, timings: List[float]) -> str: + total = sum(timings) + mean = total / len(timings) + return ( + f"{label:<28} total {total:6.3f} s | mean {mean * 1000:7.2f} ms | " + f"min {min(timings) * 1000:7.2f} ms | max {max(timings) * 1000:7.2f} ms" + ) + + +def main() -> None: + project_dir = Path(__file__).resolve().parents[2] + data_dir = project_dir / "data" / "TUDELFT_V3_KITE" + ml_models_dir = project_dir / "data" / "ml_models" + + configs = [ + ( + "Breukels regression", + data_dir + / "CAD_derived_geometry" + / "aero_geometry_CAD_breukels_regression.yaml", + None, + ), + ( + "Masure ML regression", + data_dir + / "CAD_derived_geometry" + / "aero_geometry_CAD_masure_regression.yaml", + ml_models_dir, + ), + ] + + n_panels = 25 + spanwise_distribution = "uniform" + Umag = 10.0 + side_slip = 0.0 + yaw_rate = 0.0 + alpha_values = [0.0, 4.0, 8.0, 12.0, 14.0] + + print("TUDELFT_V3_KITE benchmark") + print( + f"Settings -> n_panels={n_panels}, Umag={Umag} m/s, " + f"alpha sweep={alpha_values}, side_slip={side_slip}, yaw_rate={yaw_rate}" + ) + + for label, config_path, model_dir in configs: + print(f"\n=== {label} ===") + + baseline_body = instantiate_body( + config_path, + n_panels, + spanwise_distribution, + model_dir, + ) + baseline_solver = Solver(gamma_initial_distribution_type="elliptical") + baseline_times, baseline_results = run_sweep( + baseline_body, + baseline_solver, + Umag, + alpha_values, + side_slip, + yaw_rate, + warm_start=False, + ) + + iterative_body = instantiate_body( + config_path, + n_panels, + spanwise_distribution, + model_dir, + ) + iterative_solver = Solver(gamma_initial_distribution_type="previous") + iterative_times, iterative_results = run_sweep( + iterative_body, + iterative_solver, + Umag, + alpha_values, + side_slip, + yaw_rate, + warm_start=True, + ) + + print(summarise_timings("Baseline (elliptical)", baseline_times)) + print(summarise_timings("Iterative (previous)", iterative_times)) + + total_baseline = sum(baseline_times) + total_iterative = sum(iterative_times) + if total_baseline > 0: + delta = total_baseline - total_iterative + pct = delta / total_baseline * 100.0 + print(f"Speed-up from warm start: {delta:6.3f} s ({pct:4.1f}%)") + + print("\n alpha [deg] | baseline [ms] | iterative [ms] | Δ [ms]") + for alpha, base_t, iter_t in zip(alpha_values, baseline_times, iterative_times): + delta_ms = (base_t - iter_t) * 1000.0 + print( + f"{alpha:10.1f} | {base_t * 1000:13.2f} | " + f"{iter_t * 1000:13.2f} | {delta_ms:7.2f}" + ) + + final_iter = iterative_results[-1] + print( + f"Final coefficients at α={final_iter['alpha']:.1f}° -> " + f"CL={final_iter['cl']:.3f}, CD={final_iter['cd']:.3f}, CS={final_iter['cs']:.3f}" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/TUDELFT_V3_KITE/compare_against_literature.py b/examples/TUDELFT_V3_KITE/compare_against_literature.py new file mode 100644 index 0000000..e72e3ea --- /dev/null +++ b/examples/TUDELFT_V3_KITE/compare_against_literature.py @@ -0,0 +1,216 @@ +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 + """ + ### defining paths + PROJECT_DIR = Path(__file__).resolve().parents[2] + + ### defining settings + n_panels = 50 + spanwise_panel_distribution = "uniform" + 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 = ( + Path(PROJECT_DIR) / "data" / "TUDELFT_V3_KITE" / "CAD_derived_geometry" + ) + body_aero_CAD_CFD_polars = BodyAerodynamics.instantiate( + n_panels=n_panels, + 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 / "aero_geometry_CAD_CFD_polars.yaml"), + spanwise_panel_distribution=spanwise_panel_distribution, + 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 / "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 / "aero_geometry_CAD_masure_regression.yaml" + ), + ml_models_dir=(Path(PROJECT_DIR) / "data" / "ml_models"), + spanwise_panel_distribution=spanwise_panel_distribution, + ) + + # 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_CAD_CFD_polars.va_initialize(Umag, angle_of_attack, side_slip, yaw_rate) + body_aero_CAD_CFD_polars_with_bridles.va_initialize( + Umag, angle_of_attack, side_slip, yaw_rate + ) + body_aero_CAD_neuralfoil.va_initialize(Umag, angle_of_attack, side_slip, yaw_rate) + body_aero_masure_regression.va_initialize( + Umag, angle_of_attack, side_slip, yaw_rate + ) + + # 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_alpha = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "CFD_RANS_Rey_10e5_Poland2025_alpha_sweep_beta_0.csv" + ) + path_windtunnel_alpha = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "WindTunnel_Re5e5_alpha_sweep_beta_0_Poland2025.csv" + ) + plot_polars( + solver_list=[ + solver_base_version, + solver_base_version, + solver_base_version, + solver_base_version, + ], + body_aero_list=[ + body_aero_CAD_CFD_polars, + body_aero_CAD_CFD_polars_with_bridles, + body_aero_CAD_neuralfoil, + body_aero_masure_regression, + ], + label_list=[ + "VSM CAD CFD Polars", + "VSM CAD CFD Polars with Bridles", + "VSM CAD NeuralFoil", + "VSM CAD Masure Regression", + path_cfd_alpha.stem, + path_windtunnel_alpha.stem, + ], + literature_path_list=[path_cfd_alpha, path_windtunnel_alpha], + angle_range=[0, 5, 8, 10, 12, 15, 20, 25], + 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, + ) + + # Step 5b: Plot beta sweep (side slip) + path_cfd_beta = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "CFD_RANS_Re1e6_beta_sweep_alpha_13_Vire2022_CorrectedByPoland2025.csv" + ) + path_windtunnel_beta = ( + Path(PROJECT_DIR) + / "data" + / "TUDELFT_V3_KITE" + / "3D_polars_literature" + / "WindTunnel_Re5e5_beta_sweep_alpha_13_Poland2025.csv" + ) + plot_polars( + solver_list=[ + solver_base_version, + solver_base_version, + solver_base_version, + solver_base_version, + ], + body_aero_list=[ + body_aero_CAD_CFD_polars, + body_aero_CAD_CFD_polars_with_bridles, + body_aero_CAD_neuralfoil, + body_aero_masure_regression, + ], + label_list=[ + "VSM CAD CFD Polars", + "VSM CAD CFD Polars with Bridles", + "VSM CAD NeuralFoil", + "VSM CAD Masure Regression", + path_cfd_beta.stem, + path_windtunnel_beta.stem, + ], + literature_path_list=[path_cfd_beta, path_windtunnel_beta], + angle_range=[0, 3, 6, 9, 12], + angle_type="side_slip", + angle_of_attack=6.8, + side_slip=0, + yaw_rate=0, + Umag=3.15, + title="betasweep", + data_type=".pdf", + save_path=Path(save_folder), + is_save=True, + is_show=True, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py b/examples/TUDELFT_V3_KITE/fit_polynomials_polars.py index fce951d..d7f9501 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" - / "aero_geometry_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/examples/TUDELFT_V3_KITE/plotting_the_geometry.py b/examples/TUDELFT_V3_KITE/plotting_the_geometry.py new file mode 100644 index 0000000..ed13fbd --- /dev/null +++ b/examples/TUDELFT_V3_KITE/plotting_the_geometry.py @@ -0,0 +1,143 @@ +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])) + + # 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_CAD_CFD_polars = BodyAerodynamics.instantiate( + n_panels=n_panels, + 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 / "aero_geometry_CAD_CFD_polars.yaml"), + spanwise_panel_distribution=spanwise_panel_distribution, + 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 / "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 / "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_CAD_CFD_polars.va_initialize(Umag, angle_of_attack, side_slip, yaw_rate) + body_aero_CAD_CFD_polars_with_bridles.va_initialize( + Umag, angle_of_attack, side_slip, yaw_rate + ) + body_aero_CAD_neuralfoil.va_initialize(Umag, angle_of_attack, side_slip, yaw_rate) + body_aero_masure_regression.va_initialize( + Umag, angle_of_attack, side_slip, yaw_rate + ) + + # Step 3: Plot the kite geometry using Matplotlib + """ + Visualize the panel mesh, control points, and aerodynamic centers for the selected BodyAerodynamics object. + """ + plot_geometry( + body_aero_CAD_CFD_polars, + title="TUDELFT_V3_KITE", + data_type=".pdf", + save_path=".", + is_save=False, + is_show=True, + ) + + # Step 4: Create an interactive plot using Plotly + """ + Generate an interactive 3D plot for the selected BodyAerodynamics object. + This allows for interactive exploration of the geometry and panel arrangement. + """ + interactive_plot( + body_aero_CAD_CFD_polars_with_bridles, + vel=Umag, + angle_of_attack=angle_of_attack, + side_slip=side_slip, + yaw_rate=yaw_rate, + is_with_aerodynamic_details=True, + title="TUDELFT_V3_KITE", + is_with_bridles=True, + ### uncomment the lines below to save + # is_save=True, + # save_path=Path(PROJECT_DIR) + # / "results" + # / "TUDELFT_V3_KITE" + # / "interactive_plot.html", + ) + + +if __name__ == "__main__": + main() diff --git a/examples/TUDELFT_V3_KITE/tutorial.py b/examples/TUDELFT_V3_KITE/tutorial.py index 783d65f..992e452 100644 --- a/examples/TUDELFT_V3_KITE/tutorial.py +++ b/examples/TUDELFT_V3_KITE/tutorial.py @@ -64,12 +64,12 @@ def main(): ) body_aero_CAD_CFD_polars_with_bridles = BodyAerodynamics.instantiate( n_panels=n_panels, - file_path=(cad_derived_geometry_dir / "aero_geometry_CAD_CFD_polars.yaml"), + # file_path=(cad_derived_geometry_dir / "config_kite_CAD_CFD_polars.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_manually_adjusted.yaml" ), - ml_models_dir=(Path(PROJECT_DIR) / "data" / "ml_models"), ) body_aero_CAD_neuralfoil = BodyAerodynamics.instantiate( n_panels=n_panels, diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index e8051f6..5099d57 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -839,6 +839,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 = ( @@ -855,7 +856,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, @@ -867,11 +868,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 e2ef621..5a1459e 100644 --- a/src/VSM/plotting.py +++ b/src/VSM/plotting.py @@ -560,6 +560,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 index 8fbe1e4..a4f5dac 100644 --- a/src/VSM/stability_derivatives.py +++ b/src/VSM/stability_derivatives.py @@ -374,3 +374,196 @@ def map_derivatives_to_aircraft_frame( 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