From a0311dfb92a1f7c6f4bc3b8ee550ccfce3ab4adf Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Thu, 6 Mar 2025 17:08:56 +0100 Subject: [PATCH 1/2] add initial test case --- test/verification_cases/v3/CFD_data.csv | 14 ++ .../v3/RANS_CD_alpha_struts.csv | 14 ++ .../v3/RANS_CL_alpha_struts.csv | 13 + test/verification_cases/v3/test_v3.jl | 101 ++++++++ test/verification_cases/v3/test_v3.py | 236 ++++++++++++++++++ 5 files changed, 378 insertions(+) create mode 100644 test/verification_cases/v3/CFD_data.csv create mode 100644 test/verification_cases/v3/RANS_CD_alpha_struts.csv create mode 100644 test/verification_cases/v3/RANS_CL_alpha_struts.csv create mode 100644 test/verification_cases/v3/test_v3.jl create mode 100644 test/verification_cases/v3/test_v3.py diff --git a/test/verification_cases/v3/CFD_data.csv b/test/verification_cases/v3/CFD_data.csv new file mode 100644 index 00000000..76e77152 --- /dev/null +++ b/test/verification_cases/v3/CFD_data.csv @@ -0,0 +1,14 @@ +CL,CD,aoa +-0.09999989561132601,0.07311065309660152,-4.899328101310155 +0.08648652661186279,0.052945258966302186,-1.9798614727392547 +0.46486494135386314,0.04843192628324135,3.053694290704481 +0.6756755947979644,0.06421637996302657,6.04026783780421 +0.8621622026010177,0.08170237956021868,9.02684599326258 +1.0324324706769314,0.10600706236966365,11.979866081097894 +1.13243245907819,0.12310883205338634,13.95973216219579 +1.237837878590173,0.14413199452399744,16.006712074360475 +1.2864865266118626,0.15490563256819245,17.01342184454163 +1.1081080886723789,0.19281212679574952,18.02013161472279 +1.13243245907819,0.21441074900333865,19.999997695820685 +1.1270270279674657,0.21481480068630648,21.946305709295526 +1.1459459672625518,0.27283952501085074,24.060399452526994 diff --git a/test/verification_cases/v3/RANS_CD_alpha_struts.csv b/test/verification_cases/v3/RANS_CD_alpha_struts.csv new file mode 100644 index 00000000..b2e2b766 --- /dev/null +++ b/test/verification_cases/v3/RANS_CD_alpha_struts.csv @@ -0,0 +1,14 @@ +-5.037036895751953, 0.07407404581705739 +-2.0370356241861995, 0.05308642917209197 +-0.0370356241861995, 0.048148176405165076 +3, 0.048148176405165076 +6.037038167317707, 0.06419754028320312 +9.000000000000004, 0.08148142496744797 +12, 0.10617277357313382 +14.000000000000005, 0.12345682779947921 +16.037038167317714, 0.14444444444444454 +17.074073791503913, 0.15555547078450518 +18.000002543131515, 0.19259253607855903 +20.037035624186203, 0.21481480068630648 +22.037035624186203, 0.21481480068630648 +24.037040710449222, 0.27283952501085074 \ No newline at end of file diff --git a/test/verification_cases/v3/RANS_CL_alpha_struts.csv b/test/verification_cases/v3/RANS_CL_alpha_struts.csv new file mode 100644 index 00000000..c639ef3c --- /dev/null +++ b/test/verification_cases/v3/RANS_CL_alpha_struts.csv @@ -0,0 +1,13 @@ +-4.899328101310155, -0.09999989561132601 +-1.9798614727392547, 0.08648652661186279 +3.053694290704481, 0.46486494135386314 +6.04026783780421, 0.6756755947979644 +9.02684599326258, 0.8621622026010177 +11.979866081097894, 1.0324324706769314 +13.95973216219579, 1.13243245907819 +16.006712074360475, 1.237837878590173 +17.01342184454163, 1.2864865266118626 +18.02013161472279, 1.1081080886723789 +19.999997695820685, 1.13243245907819 +21.946305709295526, 1.1270270279674657 +24.060399452526994, 1.1459459672625518 \ No newline at end of file diff --git a/test/verification_cases/v3/test_v3.jl b/test/verification_cases/v3/test_v3.jl new file mode 100644 index 00000000..be99080e --- /dev/null +++ b/test/verification_cases/v3/test_v3.jl @@ -0,0 +1,101 @@ +function get_CAD_matching_uri() + CAD_matching_pos_values = + [0.000000e00 0.000000e00 0.000000e00 + 1.538773e03 4.113307e03 5.530496e03 + -1.762300e01 3.967978e03 6.471622e03 + -2.376830e02 3.134335e03 7.476759e03 + -3.837330e02 1.959729e03 8.078914e03 + -4.565620e02 6.642520e02 8.339101e03 + -4.565620e02 -6.642520e02 8.339101e03 + -3.837330e02 -1.959729e03 8.078914e03 + -2.376830e02 -3.134335e03 7.476759e03 + -1.762300e01 -3.967978e03 6.471622e03 + 1.538773e03 -4.113307e03 5.530496e03 + 1.703467e03 -3.955506e03 6.467819e03 + 2.002516e03 -3.116753e03 7.454254e03 + 2.110145e03 -1.946587e03 8.041739e03 + 2.158559e03 -6.607210e02 8.294064e03 + 2.158559e03 6.607210e02 8.294064e03 + 2.110145e03 1.946587e03 8.041739e03 + 2.002516e03 3.116753e03 7.454254e03 + 1.703467e03 3.955506e03 6.467819e03 + 8.595800e02 4.139660e03 5.654227e03 + 8.595800e02 -4.139660e03 5.654227e03 + 3.138480e02 0.000000e00 1.252129e03 + 3.875460e02 6.909170e02 1.479174e03 + 1.053321e03 1.772499e03 4.343864e03 + 1.441205e03 2.708913e03 5.601785e03 + 1.528031e03 1.338349e03 5.966178e03 + 3.875460e02 -6.909170e02 1.479174e03 + 1.053321e03 -1.772499e03 4.343864e03 + 1.441205e03 -2.708913e03 5.601785e03 + 1.528031e03 -1.338349e03 5.966178e03 + -6.552600e01 1.321471e03 4.213046e03 + -4.248900e01 2.046976e03 5.525725e03 + -9.173800e01 1.262274e03 5.961848e03 + -6.552600e01 -1.321471e03 4.213046e03 + -4.248900e01 -2.046976e03 5.525725e03 + -9.173800e01 -1.262274e03 5.961848e03] + return CAD_matching_pos_values +end + +function get_v3_case_params() + wing_type = "LEI_kite" + dist = "lin" + N_split = 4 + aoas = collect(-4:2:22) + Umag = 22 + # convergence criteria + max_iterations = 1500 + allowed_error = 1e-5 + relaxation_factor = 0.03 + core_radius_fraction = 1e-20 + + # Wing geometry + coord_struc = get_CAD_matching_uri() + # coord = thesis_functions.struct2aero_geometry(coord_struc) / 1000 + + # N = len(coord) // 2 + + # # LE thickness at each section [m] + # # 10 sections + # LE_thicc = 0.1 + + # # Camber for each section (ct in my case) + # camber = 0.095 + + # # Refine structrural mesh into more panels + # coord = thesis_functions.refine_LEI_mesh(coord, N - 1, N_split) + # N = int(len(coord) / 2) # Number of section after refining the mesh + + # # Definition of airfoil coefficients + # # Based on Breukels (2011) correlation model + # aoas_for_polar = np.arange(-80, 80, 0.1) + # data_airf = np.empty((len(aoas_for_polar), 4)) + # for j in range(len(aoas_for_polar)) + # alpha = aoas_for_polar[j] + # Cl, Cd, Cm = thesis_functions.LEI_airf_coeff(LE_thicc, camber, alpha) + # data_airf[j, 0] = alpha + # data_airf[j, 1] = Cl + # data_airf[j, 2] = Cd + # data_airf[j, 3] = Cm + # end + + # Atot = test_utils.calculate_projected_area(coord) + # coord_input_params = [coord, LE_thicc, camber] + # case_parameters = [ + # coord_input_params, + # aoas, + # wing_type, + # Umag, + # 0, + # Atot, + # max_iterations, + # allowed_error, + # relaxation_factor, + # core_radius_fraction, + # data_airf, + # ] + # + # return case_parameters +end \ No newline at end of file diff --git a/test/verification_cases/v3/test_v3.py b/test/verification_cases/v3/test_v3.py new file mode 100644 index 00000000..f26148f4 --- /dev/null +++ b/test/verification_cases/v3/test_v3.py @@ -0,0 +1,236 @@ +# -*- coding: utf-8 -*- +import numpy as np +import logging +import sys +import os +import matplotlib.pyplot as plt +import pytest + +# Go back to root folder +root_path = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +sys.path.insert(0, root_path) +import tests.utils as test_utils +import tests.thesis_functions_oriol_cayon as thesis_functions + + +@pytest.fixture(autouse=True) +def change_test_dir(request): + # Get the directory of the current test file + test_dir = os.path.dirname(request.fspath) + # Change the working directory to the test file's directory + os.chdir(test_dir) + + +def get_v3_case_params(): + + wing_type = "LEI_kite" + dist = "lin" + N_split = 4 + aoas = np.arange(-4, 24, 2) + Umag = 22 + # convergence criteria + max_iterations = 1500 + allowed_error = 1e-5 + relaxation_factor = 0.03 + core_radius_fraction = 1e-20 + + # Wing geometry + coord_struc = thesis_functions.get_CAD_matching_uri() + coord = thesis_functions.struct2aero_geometry(coord_struc) / 1000 + + N = len(coord) // 2 + + # LE thickness at each section [m] + # 10 sections + LE_thicc = 0.1 + + # Camber for each section (ct in my case) + camber = 0.095 + + # Refine structrural mesh into more panels + coord = thesis_functions.refine_LEI_mesh(coord, N - 1, N_split) + N = int(len(coord) / 2) # Number of section after refining the mesh + + # Definition of airfoil coefficients + # Based on Breukels (2011) correlation model + aoas_for_polar = np.arange(-80, 80, 0.1) + data_airf = np.empty((len(aoas_for_polar), 4)) + for j in range(len(aoas_for_polar)): + alpha = aoas_for_polar[j] + Cl, Cd, Cm = thesis_functions.LEI_airf_coeff(LE_thicc, camber, alpha) + data_airf[j, 0] = alpha + data_airf[j, 1] = Cl + data_airf[j, 2] = Cd + data_airf[j, 3] = Cm + + Atot = test_utils.calculate_projected_area(coord) + coord_input_params = [coord, LE_thicc, camber] + case_parameters = [ + coord_input_params, + aoas, + wing_type, + Umag, + 0, + Atot, + max_iterations, + allowed_error, + relaxation_factor, + core_radius_fraction, + data_airf, + ] + + return case_parameters + + +def test_v3(): + case_params = get_v3_case_params() + # making sure not too many points are tested + case_params[1] = np.deg2rad(np.array([4, 12])) + # changing wing_type to take the polars and not polynomial directly + case_params[2] = "LEI_kite_polars" + # comparison solution + aoas = case_params[1] + + ### COMPARING FROM POLARS + # OLD numerical + CL_LLT, CD_LLT, CL_VSM, CD_VSM, gamma_LLT, gamma_VSM = ( + test_utils.calculate_old_for_alpha_range(case_params) + ) + # NEW numerical + ( + CL_LLT_new, + CD_LLT_new, + CL_VSM_new, + CD_VSM_new, + gamma_LLT_new, + gamma_VSM_new, + panel_y, + AR, + ) = test_utils.calculate_new_for_alpha_range( + case_params, + is_plotting=False, + ) + # checking LTT old close to LLT new + for gamma_LLT_i, gamma_LLT_new_i in zip(gamma_LLT, gamma_LLT_new): + assert np.allclose(gamma_LLT_i, gamma_LLT_new_i, atol=1e-2) + assert np.allclose(CL_LLT, CL_LLT_new, atol=1e-3) + assert np.allclose(CD_LLT, CD_LLT_new, atol=1e-4) + + # checking VSMs to be close to one another + assert np.allclose(CL_VSM, CL_VSM_new, atol=1e-3) + assert np.allclose(CD_VSM, CD_VSM_new, atol=1e-4) + + ################################################## + ### COMPARING FROM POLYNOMIAL + case_params = get_v3_case_params() + # making sure not too many points are tested + case_params[1] = np.deg2rad(np.array([3, 6, 9])) + # changing wing_type to take the polars and not polynomial directly + case_params[2] = "LEI_kite" + # comparison solution + aoas = case_params[1] + + # changing wing_type to take the polars and not polynomial directly + + # OLD numerical + CL_LLT, CD_LLT, CL_VSM, CD_VSM, gamma_LLT, gamma_VSM = ( + test_utils.calculate_old_for_alpha_range(case_params) + ) + # NEW numerical + ( + CL_LLT_new, + CD_LLT_new, + CL_VSM_new, + CD_VSM_new, + gamma_LLT_new, + gamma_VSM_new, + panel_y, + AR, + ) = test_utils.calculate_new_for_alpha_range( + case_params, + is_plotting=False, + ) + + # checking LTT old close to LLT new + assert np.allclose(CL_LLT, CL_LLT_new, atol=2e-2) + assert np.allclose(CD_LLT, CD_LLT_new, atol=2e-3) + + # checking VSMs to be close to one another + for gamma_VSM_i, gamma_VSM_new_i in zip(gamma_VSM, gamma_VSM_new): + assert np.allclose(gamma_VSM_i, gamma_VSM_new_i, atol=1e-2) + assert np.allclose(CL_VSM, CL_VSM_new, atol=2e-2) + assert np.allclose(CD_VSM, CD_VSM_new, atol=2e-3) + + # comparing solution + CL_struts = np.loadtxt("./CFD_data/RANS_CL_alpha_struts.csv", delimiter=",") + CD_struts = np.loadtxt("./CFD_data/RANS_CD_alpha_struts.csv", delimiter=",") + + CL_CFD = CL_struts[:, 1] + alpha_CFD = CL_struts[:, 0] + CD_CFD = np.interp(alpha_CFD, CD_struts[:, 0], CD_struts[:, 1]) + + CL_CFD_at_new_alphas = [] + CD_CFD_at_new_alphas = [] + for alpha in aoas: + alpha = np.rad2deg(alpha) + CL_CFD_at_new_alphas.append(np.interp(alpha, alpha_CFD, CL_CFD)) + CD_CFD_at_new_alphas.append(np.interp(alpha, alpha_CFD, CD_CFD)) + + # checking new VSM close to Maneia + assert np.allclose(CL_CFD_at_new_alphas, CL_VSM_new, atol=1e-1) + assert np.allclose(CD_CFD_at_new_alphas, CD_VSM_new, atol=1e-2) + + +if __name__ == "__main__": + + case_params = get_v3_case_params() + + aoas = np.deg2rad(np.linspace(0, 25, 30)) + aoas = np.deg2rad(np.array([30, 32])) + case_params[1] = aoas + # comparing solution + CL_struts = np.loadtxt("./CFD_data/RANS_CL_alpha_struts.csv", delimiter=",") + CD_struts = np.loadtxt("./CFD_data/RANS_CD_alpha_struts.csv", delimiter=",") + + CL_CFD = CL_struts[:, 1] + aoas_CFD = CL_struts[:, 0] + CD_CFD = np.interp(aoas_CFD, CD_struts[:, 0], CD_struts[:, 1]) + polars_CFD = np.vstack((aoas_CFD, CL_CFD, CD_CFD)).T + # OLD numerical + CL_LLT, CD_LLT, CL_VSM, CD_VSM, gamma_LLT, gamma_VSM = ( + test_utils.calculate_old_for_alpha_range(case_params) + ) + # NEW numerical + ( + CL_LLT_new, + CD_LLT_new, + CL_VSM_new, + CD_VSM_new, + gamma_LLT_new, + gamma_VSM_new, + panel_y, + AR, + ) = test_utils.calculate_new_for_alpha_range( + case_params, + is_plotting=False, + ) + + test_utils.plotting_CL_CD_gamma_LLT_VSM_old_new_comparison( + panel_y=panel_y, + AR=AR, + wing_type="LEI_kite", + aoas=[polars_CFD[:, 0], aoas], + CL_list=[polars_CFD[:, 1], CL_LLT, CL_LLT_new, CL_VSM, CL_VSM_new], + CD_list=[polars_CFD[:, 2], CD_LLT, CD_LLT_new, CD_VSM, CD_VSM_new], + gamma_list=[gamma_LLT, gamma_LLT_new, gamma_VSM, gamma_VSM_new], + labels=["RANS [Lebesque]", "LLT", "LLT_new", "VSM", "VSM_new"], + ) + labels = ["Polars CFD", "LLT", "LLT_new", "VSM", "VSM_new"] + CL_list = [polars_CFD[:, 1], CL_LLT, CL_LLT_new, CL_VSM, CL_VSM_new] + CD_list = [polars_CFD[:, 2], CD_LLT, CD_LLT_new, CD_VSM, CD_VSM_new] + for i, aoa in enumerate(aoas): + print(f"aoa = {np.rad2deg(aoa)}") + for label, CD, CL in zip(labels, CD_list, CL_list): + print(f"{label}: CL = {CL}, CD = {CD}") + + plt.show() From 4e1607490f1a7a155de334feab7c03b3349f119b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Thu, 6 Mar 2025 17:22:12 +0100 Subject: [PATCH 2/2] some progress --- test/verification_cases/v3/test_v3.jl | 41 +++++++++++++++++++++------ 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/test/verification_cases/v3/test_v3.jl b/test/verification_cases/v3/test_v3.jl index be99080e..6566f9e8 100644 --- a/test/verification_cases/v3/test_v3.jl +++ b/test/verification_cases/v3/test_v3.jl @@ -38,6 +38,31 @@ function get_CAD_matching_uri() -9.173800e01 -1.262274e03 5.961848e03] return CAD_matching_pos_values end +function struct2aero_geometry(coord_struc) + coord = zeros(20, 3) + coord[1, :] = coord_struc[21, :] + coord[2, :] = coord_struc[11, :] + coord[3, :] = coord_struc[10, :] + coord[4, :] = coord_struc[12, :] + coord[5, :] = coord_struc[9, :] + coord[6, :] = coord_struc[13, :] + coord[7, :] = coord_struc[8, :] + coord[8, :] = coord_struc[14, :] + coord[9, :] = coord_struc[7, :] + coord[10, :] = coord_struc[15, :] + coord[11, :] = coord_struc[6, :] + coord[12, :] = coord_struc[16, :] + coord[13, :] = coord_struc[5, :] + coord[14, :] = coord_struc[17, :] + coord[15, :] = coord_struc[4, :] + coord[16, :] = coord_struc[18, :] + coord[17, :] = coord_struc[3, :] + coord[18, :] = coord_struc[19, :] + coord[19, :] = coord_struc[20, :] + coord[20, :] = coord_struc[2, :] + return coord +end + function get_v3_case_params() wing_type = "LEI_kite" @@ -53,18 +78,18 @@ function get_v3_case_params() # Wing geometry coord_struc = get_CAD_matching_uri() - # coord = thesis_functions.struct2aero_geometry(coord_struc) / 1000 + coord = struct2aero_geometry(coord_struc) / 1000 - # N = len(coord) // 2 + N = length(coord) รท 2 - # # LE thickness at each section [m] - # # 10 sections - # LE_thicc = 0.1 + # LE thickness at each section [m] + # 10 sections + LE_thicc = 0.1 - # # Camber for each section (ct in my case) - # camber = 0.095 + # Camber for each section (ct in my case) + camber = 0.095 - # # Refine structrural mesh into more panels + # # Refine structural mesh into more panels # coord = thesis_functions.refine_LEI_mesh(coord, N - 1, N_split) # N = int(len(coord) / 2) # Number of section after refining the mesh