From b784646a50e55178153dd106c3434dd93591d92f Mon Sep 17 00:00:00 2001 From: Samuel Ayala Date: Tue, 21 May 2024 14:50:50 -0400 Subject: [PATCH] better theodorsen plotting --- data/theodorsen.py | 102 ++++++++++++++++----------------- tests/test_uvlm_theodorsen.cpp | 68 +++++++++++----------- 2 files changed, 85 insertions(+), 85 deletions(-) diff --git a/data/theodorsen.py b/data/theodorsen.py index 9292ace..7e93d7b 100644 --- a/data/theodorsen.py +++ b/data/theodorsen.py @@ -21,6 +21,24 @@ def theo_fun(k): C = H1 / (H1 + 1.0j * H0) return C +def uvlm_data(filename): + uvlm_cl = [] + uvlm_t = [] + uvlm_z = [] + uvlm_alpha = [] + with open(f"build/windows/x64/release/{filename}.txt", "r") as f: + k = float(f.readline()) # get reduced frequency of the simulation + for line in f: + time_step, z, cl, alpha = line.split() + uvlm_t.append(float(time_step)) + uvlm_z.append(float(z)) + uvlm_cl.append(float(cl)) + uvlm_alpha.append(float(alpha)) + + uvlm_alpha = np.array(uvlm_alpha) + uvlm_cl = np.array(uvlm_cl) + return k, uvlm_t, uvlm_z, uvlm_cl, uvlm_alpha + # Some info in Katz Plotkin p414 (eq 13.73a) # Jone's approximation of Wagner function b0 = 1 @@ -38,33 +56,10 @@ def theo_fun(k): a = ar / c # full span pitch_axis = -0.5 # leading edge +# Theodorsen numerical simulation param cycles = 3 # number of periods - -uvlm_cl = [] -uvlm_t = [] -uvlm_z = [] -uvlm_alpha = [] -with open("build/windows/x64/release/cl_data.txt", "r") as f: - freq = float(f.readline()) # get reduced frequency of the simulation - for line in f: - time_step, z, cl, alpha = line.split() - uvlm_t.append(float(time_step)) - uvlm_z.append(float(z)) - uvlm_cl.append(float(cl)) - uvlm_alpha.append(float(alpha)) - -n = len(uvlm_t) # number of time steps -uvlm_cycle_idx = int((1 - 1 / cycles) * n - 1) -uvlm_alpha = np.array(uvlm_alpha) -uvlm_cl = np.array(uvlm_cl) - -amplitudes = [0.1] -reduced_frequencies = [freq] - -t_final = cycles * 2 * np.pi / (reduced_frequencies[0] * 2.0 * u_inf / c) # 4 cycles -nb_pts = 500 +nb_pts = 1000 cycle_idx = int((1 - 1 / cycles) * nb_pts - 1) -vec_t = np.linspace(0, t_final, nb_pts) fig, axs = plt.subplot_mosaic( [["time"], ["heave"]], # Disposition des graphiques @@ -72,22 +67,34 @@ def theo_fun(k): figsize=(11, 6), # Ajuster la taille de la figure (x,y) ) -for amp, k in zip(amplitudes, reduced_frequencies): - # heave parameters - amplitude = amp - omega = k * 2.0 * u_inf / c # pitch frequency +files = [ + "cl_data_01", + "cl_data_03", + "cl_data_05", + "cl_data_07", +] + +for file in files: + k, uvlm_t, uvlm_z, uvlm_cl, uvlm_alpha = uvlm_data(file) + t_final = cycles * 2 * np.pi / (k * 2.0 * u_inf / c) + omega = k * 2.0 * u_inf / c # frequency + + vec_t = np.linspace(0, t_final, nb_pts) + + n = len(uvlm_t) # number of time steps + uvlm_cycle_idx = int((1 - 1 / cycles) * n - 1) # sudden acceleration # def pitch(t): return np.radians(5) # def heave(t): return 0 # pure heaving - # def pitch(t): return 0 - # def heave(t): return -amplitude * np.sin(omega * t) + def pitch(t): return 0 + def heave(t): return -0.1 * np.sin(omega * t) # pure pitching - def pitch(t): return np.radians(np.sin(omega * t)) - def heave(t): return 0 + # def pitch(t): return np.radians(np.sin(omega * t)) + # def heave(t): return 0 def w(s: float): return u_inf * pitch(s) + derivative(heave, s) + b * (0.5 - pitch_axis) * derivative(pitch, s) @@ -107,30 +114,19 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation return (L_m + L_c) / (0.5 * rho * u_inf * u_inf * c) cl_theo = np.array([cl_theodorsen(ti) for ti in vec_t]) - coord_z = np.array([-heave(ti) / c for ti in vec_t]) + h = np.array([-heave(ti) / c for ti in vec_t]) angle = np.array([np.degrees(pitch(ti)) for ti in vec_t]) - axs["time"].plot(vec_t, cl_theo, label=f"k={k}") - axs["heave"].plot(angle[cycle_idx:], cl_theo[cycle_idx:], label=f"k={k}") - -analytical_cl = np.array([np.interp(ut, vec_t, cl_theo) for ut in uvlm_t[uvlm_cycle_idx:]]) -difference = uvlm_cl[uvlm_cycle_idx:] - analytical_cl -error = np.sqrt(np.dot(difference, difference) / (n-uvlm_cycle_idx)) -print(f"Error: {100 * error / np.max(np.abs(analytical_cl)):.3f}%", ) - -# katz_t = [] -# katz_cl = [] -# with open("data/katz/katz_cl.txt", "r") as f: -# for line in f: -# t, cl = line.split() -# katz_t.append(float(t)) -# katz_cl.append(float(cl)) + plotc, = axs["time"].plot(vec_t, cl_theo, label=f"Theodorsen (k={k})") + axs["heave"].plot(h[cycle_idx:], cl_theo[cycle_idx:], label=f"Theodorsen (k={k})") -# axs["time"].scatter(katz_t[1:], katz_cl[1:], label="Katz", facecolors='none', edgecolors='g', s=20) + axs["time"].scatter(uvlm_t, uvlm_cl, label=f"UVLM (k={k})", facecolors='none', edgecolors=plotc.get_color(), s=20) + axs["heave"].scatter(uvlm_z[uvlm_cycle_idx:], uvlm_cl[uvlm_cycle_idx:], label=f"UVLM (k={k})", facecolors='none', edgecolors=plotc.get_color(), s=20) -axs["time"].scatter(uvlm_t, uvlm_cl, label="UVLM", facecolors='none', edgecolors='r', s=20) -# axs["time"].scatter(uvlm_t, uvlm_z, label="UVLM z", facecolors='none', edgecolors='r', s=20) -axs["heave"].scatter(uvlm_alpha[uvlm_cycle_idx:], uvlm_cl[uvlm_cycle_idx:], label="UVLM", facecolors='none', edgecolors='r', s=20) + analytical_cl = np.array([np.interp(ut, vec_t, cl_theo) for ut in uvlm_t[uvlm_cycle_idx:]]) + difference = uvlm_cl[uvlm_cycle_idx:] - analytical_cl + error = np.sqrt(np.dot(difference, difference) / (n-uvlm_cycle_idx)) + print(f"Freq: {k}, Error: {100 * error / np.max(np.abs(analytical_cl)):.3f}%", ) # axs["time"].plot(vec_t, [0.548311] * len(vec_t), label="VLM (alpha=5)") diff --git a/tests/test_uvlm_theodorsen.cpp b/tests/test_uvlm_theodorsen.cpp index 402856e..d0c49b9 100644 --- a/tests/test_uvlm_theodorsen.cpp +++ b/tests/test_uvlm_theodorsen.cpp @@ -95,7 +95,7 @@ int main() { const f32 cycles = 3.0f; const f32 u_inf = 1.0f; // freestream velocity const f32 amplitude = 0.1f; // amplitude of the wing motion - const f32 k = 0.2; // reduced frequency + const f32 k = 0.7; // reduced frequency const f32 omega = k * 2.0f * u_inf / (2*b); const f32 t_final = cycles * 2.0f * PI_f / omega; // 4 periods //const f32 t_final = 5.0f; @@ -111,36 +111,36 @@ int main() { ); // Periodic heaving - // kinematics.add([=](f32 t) { - // return linalg::translation_matrix(linalg::alias::float3{ - // -u_inf*t, - // 0.0f, - // 0.0f - // }); - // }); - // kinematics.add([=](f32 t) { - // return linalg::translation_matrix(linalg::alias::float3{ - // 0.0f, - // 0.0f, - // amplitude * std::sin(omega * t) // heaving - // }); - // }); - - // Periodic pitching kinematics.add([=](f32 t) { return linalg::translation_matrix(linalg::alias::float3{ - -u_inf*t, // freestream + -u_inf*t, 0.0f, 0.0f }); }); kinematics.add([=](f32 t) { - return linalg::rotation_matrix( - linalg::alias::float3{0.25f, 0.0f, 0.0f}, - linalg::alias::float3{0.0f, 1.0f, 0.0f}, - to_radians(std::sin(omega * t)) - ); + return linalg::translation_matrix(linalg::alias::float3{ + 0.0f, + 0.0f, + amplitude * std::sin(omega * t) // heaving + }); }); + + // Periodic pitching + // kinematics.add([=](f32 t) { + // return linalg::translation_matrix(linalg::alias::float3{ + // -u_inf*t, // freestream + // 0.0f, + // 0.0f + // }); + // }); + // kinematics.add([=](f32 t) { + // return linalg::rotation_matrix( + // linalg::alias::float3{0.25f, 0.0f, 0.0f}, + // linalg::alias::float3{0.0f, 1.0f, 0.0f}, + // to_radians(std::sin(omega * t)) + // ); + // }); // Sudden acceleration // const f32 alpha = to_radians(5.0f); @@ -252,19 +252,23 @@ int main() { linalg::alias::float3 freestream; for (u64 idx = 0; idx < mesh->nb_panels_wing(); idx++) { const linalg::alias::float4 colloc_pt{mesh->colloc.x[idx], mesh->colloc.y[idx], mesh->colloc.z[idx], 1.0f}; - auto local_velocity = -kinematics.velocity(t, colloc_pt); - velocities.x[idx] = local_velocity.x; - velocities.y[idx] = local_velocity.y; - velocities.z[idx] = local_velocity.z; + // auto local_velocity = -kinematics.velocity(t, colloc_pt); + // velocities.x[idx] = local_velocity.x; + // velocities.y[idx] = local_velocity.y; + // velocities.z[idx] = local_velocity.z; + + velocities.x[idx] = u_inf; + velocities.y[idx] = 0.0f; + velocities.z[idx] = - omega * amplitude * std::cos(omega * t); if (idx == 0) { linalg::alias::float4 freestream_vel = -kinematics.velocity(t, colloc_pt, 1); freestream = {freestream_vel.x, freestream_vel.y, freestream_vel.z}; - std::cout << "Freestream: " << freestream << "\n"; - const f32 analytical_vel = - amplitude * omega * std::cos(omega * t); - const f32 rel_error = 100.0f * std::abs((analytical_vel - local_velocity.z) / analytical_vel); - std::cout << "vel error:" << rel_error << "%\n"; - avg_vel_error += rel_error; + // std::cout << "Freestream: " << freestream << "\n"; + // const f32 analytical_vel = - amplitude * omega * std::cos(omega * t); + // const f32 rel_error = 100.0f * std::abs((analytical_vel - local_velocity.z) / analytical_vel); + // std::cout << "vel error:" << rel_error << "%\n"; + // avg_vel_error += rel_error; } }