From 7f924562617fa3940e55eaab8f87e6ebb8f4ed68 Mon Sep 17 00:00:00 2001 From: Samuel Ayala Date: Thu, 11 Apr 2024 20:38:49 -0400 Subject: [PATCH] improvements to theodorsen we are close as the shape is good but we have some sort of negative offset with acceleration of a plate --- data/theodorsen.py | 45 +++++++++++---- mesh/infinite_rectangular_5x20.x | 98 ++++++++++++++++++++++++++++++++ tests/test_uvlm_theodorsen.cpp | 4 +- 3 files changed, 134 insertions(+), 13 deletions(-) create mode 100644 mesh/infinite_rectangular_5x20.x diff --git a/data/theodorsen.py b/data/theodorsen.py index 8047351..3ac7fe3 100644 --- a/data/theodorsen.py +++ b/data/theodorsen.py @@ -14,6 +14,9 @@ def derivative2(f, x): def solve_ivp(x0: float, s0: float, sf: float, f: callable): return spi.solve_ivp(f, [s0, sf], [x0]).y[-1] # return only the result at t=sf +def pade_approximation(k: float): + return (0.5177*k*k + 0.2752*k + 0.01576) / (k*k + 0.3414*k + 0.01582) + # Some info in Katz Plotkin p414 (eq 13.73a) # Jone's approximation of Wagner function b0 = 1 @@ -25,9 +28,13 @@ def solve_ivp(x0: float, s0: float, sf: float, f: callable): # UVLM parameters rho = 1 # fluid density u_inf = 1 # freestream -ar = 4 # aspect ratio +ar = 500 # aspect ratio b = 0.5 # half chord -a = ar / (2*b) # full span +c = 2*b # chord +a = ar / c # full span +pitch_axis = -1 # leading edge + +def atime(t: float): return 2. * u_inf * t / c amplitudes = [0.1, 0.1, 0.1] reduced_frequencies = [0.5, 0.75, 1.5] @@ -46,26 +53,42 @@ def solve_ivp(x0: float, s0: float, sf: float, f: callable): amplitude = amp / (2*b) omega = k * u_inf / (2*b) # pitch frequency - def pitch(t): return 0 - def heave(t): return amplitude * np.sin(omega * t) + # sudden acceleration + def pitch(t): return np.radians(5) + def heave(t): return 0 - def w(s: float): return u_inf * pitch(s) + derivative(heave, s) + b * (0.5 - a) * derivative(pitch, s) + # pure heaving + # def pitch(t): return 0 + # def heave(t): return amplitude * np.sin(omega * t) + + def w(s: float): + return u_inf * pitch(s) + derivative(heave, s) + b * (0.5 - pitch_axis) * derivative(pitch, s) def dx1ds(s: float, x1: float): return b1 * beta_1 * w(s) - beta_1 * x1 def dx2ds(s: float, x2: float): return b2 * beta_2 * w(s) - beta_2 * x2 - def cl_theodorsen(t: float): - L_m = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * a * derivative2(pitch, t)) - L_c = -2 * np.pi * rho * u_inf * b * ((b0 + b1 + b2) * w(t) + solve_ivp(0, 0, t, dx1ds)[-1] + solve_ivp(0, 0, t, dx2ds)[-1]) - return (L_m + L_c) / (0.5 * rho * u_inf * u_inf * a * b) - + x1_solution = spi.solve_ivp(dx1ds, [0, t_final], [0], t_eval=t) + x2_solution = spi.solve_ivp(dx2ds, [0, t_final], [0], t_eval=t) + + def x1(s: float): return np.interp(s, x1_solution.t, x1_solution.y[0]) + def x2(s: float): return np.interp(s, x2_solution.t, x2_solution.y[0]) + + def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation + L_m = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * pitch_axis * derivative2(pitch, t)) + L_c = -2 * np.pi * rho * u_inf * b * ((b0 + b1 + b2) * w(t) + x1(t) + x2(t)) + return (L_m + L_c) / (0.5 * rho * u_inf * u_inf * c) + + # def cl_theodorsen(t: float): # using Pade approximation + # return 0.5 * np.pi * (derivative2(heave, t) + derivative(pitch, t) - 0.5 * pitch_axis * derivative2(pitch, t)) + 2.0 * np.pi * (pitch(t) + derivative(heave, t) + 0.5 * derivative(pitch, t) * (0.5 - pitch_axis)) * pade_approximation(k) + # L = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * pitch_axis * derivative2(pitch, t)) + 2 * np.pi * rho * u_inf * b * (u_inf * pitch(t) + derivative(heave, t) + b * (0.5 - pitch_axis) * derivative(pitch, t)) * pade_approximation(k) + # return L / (0.5 * rho * u_inf * u_inf * a * (2*b)) + cl = np.array([cl_theodorsen(ti) for ti in t]) coord_z = np.array([heave(ti) / (2*b) for ti in t]) axs["time"].plot(t, cl, label=f"k={k}") axs["heave"].plot(coord_z[len(cl)//2:], cl[len(cl)//2:], label=f"k={k}") - uvlm_cl = [] uvlm_t = [] uvlm_z = [] diff --git a/mesh/infinite_rectangular_5x20.x b/mesh/infinite_rectangular_5x20.x new file mode 100644 index 0000000..1dde83a --- /dev/null +++ b/mesh/infinite_rectangular_5x20.x @@ -0,0 +1,98 @@ + 1 + 6 21 1 + 0.0000000000000000 0.19999999999999996 0.39999999999999991 0.59999999999999998 + 0.80000000000000004 1.0000000000000000 0.0000000000000000 0.19999999999999996 + 0.39999999999999991 0.59999999999999998 0.79999999999999993 1.0000000000000000 + 0.0000000000000000 0.19999999999999998 0.39999999999999991 0.59999999999999998 + 0.80000000000000004 1.0000000000000000 0.0000000000000000 0.19999999999999996 + 0.39999999999999991 0.60000000000000009 0.80000000000000016 1.0000000000000000 + 0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998 + 0.80000000000000004 1.0000000000000000 0.0000000000000000 0.20000000000000001 + 0.39999999999999991 0.59999999999999976 0.80000000000000004 1.0000000000000000 + 0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000001 + 0.39999999999999991 0.60000000000000020 0.79999999999999993 1.0000000000000000 + 0.0000000000000000 0.20000000000000004 0.39999999999999991 0.59999999999999998 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000007 + 0.39999999999999991 0.60000000000000020 0.80000000000000004 1.0000000000000000 + 0.0000000000000000 0.20000000000000007 0.39999999999999991 0.60000000000000009 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000007 + 0.39999999999999991 0.60000000000000009 0.79999999999999971 1.0000000000000000 + 0.0000000000000000 0.20000000000000009 0.39999999999999991 0.60000000000000009 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000009 + 0.39999999999999991 0.60000000000000009 0.79999999999999993 1.0000000000000000 + 0.0000000000000000 0.20000000000000012 0.39999999999999991 0.60000000000000009 + 0.79999999999999960 1.0000000000000000 0.0000000000000000 0.20000000000000012 + 0.39999999999999991 0.60000000000000009 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.20000000000000012 0.39999999999999991 0.60000000000000009 + 0.79999999999999982 1.0000000000000000 0.0000000000000000 0.20000000000000015 + 0.39999999999999991 0.60000000000000009 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.20000000000000015 0.39999999999999991 0.60000000000000009 + 0.79999999999999982 1.0000000000000000 0.0000000000000000 0.20000000000000018 + 0.39999999999999991 0.60000000000000009 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.20000000000000018 0.39999999999999991 0.60000000000000009 + 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 2000.0000000000000 1999.9999999999986 + 1999.9999999999973 1999.9999999999959 1999.9999999999945 1999.9999999999930 + 4000.0000000000000 3999.9999999999973 3999.9999999999945 3999.9999999999918 + 3999.9999999999891 3999.9999999999859 6000.0000000000000 6000.0000000000027 + 6000.0000000000045 6000.0000000000091 6000.0000000000118 6000.0000000000146 + 8000.0000000000000 8000.0000000000018 8000.0000000000045 8000.0000000000045 + 8000.0000000000045 8000.0000000000073 10000.000000000000 10000.000000000000 + 10000.000000000000 10000.000000000000 10000.000000000000 10000.000000000000 + 12000.000000000000 11999.999999999998 11999.999999999995 11999.999999999998 + 11999.999999999991 11999.999999999993 14000.000000000000 13999.999999999998 + 14000.000000000002 14000.000000000004 14000.000000000004 14000.000000000004 + 16000.000000000000 16000.000000000000 16000.000000000002 15999.999999999998 + 15999.999999999996 15999.999999999996 18000.000000000000 17999.999999999996 + 18000.000000000004 18000.000000000004 18000.000000000011 18000.000000000007 + 20000.000000000000 20000.000000000000 20000.000000000000 20000.000000000000 + 20000.000000000000 20000.000000000000 22000.000000000004 22000.000000000000 + 22000.000000000007 22000.000000000007 22000.000000000000 22000.000000000004 + 24000.000000000000 24000.000000000000 24000.000000000000 23999.999999999996 + 23999.999999999996 24000.000000000004 26000.000000000000 26000.000000000000 + 26000.000000000000 26000.000000000004 26000.000000000004 25999.999999999996 + 28000.000000000000 27999.999999999996 27999.999999999996 27999.999999999996 + 27999.999999999996 28000.000000000000 30000.000000000000 30000.000000000000 + 30000.000000000000 30000.000000000000 30000.000000000000 30000.000000000000 + 32000.000000000000 32000.000000000000 32000.000000000000 32000.000000000000 + 32000.000000000000 32000.000000000000 34000.000000000000 34000.000000000000 + 34000.000000000000 34000.000000000000 34000.000000000000 34000.000000000000 + 36000.000000000000 36000.000000000000 36000.000000000000 36000.000000000000 + 36000.000000000000 36000.000000000000 38000.000000000000 38000.000000000000 + 38000.000000000000 38000.000000000000 38000.000000000000 38000.000000000000 + 40000.000000000000 40000.000000000000 40000.000000000000 40000.000000000000 + 40000.000000000000 40000.000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 diff --git a/tests/test_uvlm_theodorsen.cpp b/tests/test_uvlm_theodorsen.cpp index 16f2d02..409fca4 100644 --- a/tests/test_uvlm_theodorsen.cpp +++ b/tests/test_uvlm_theodorsen.cpp @@ -64,7 +64,7 @@ void dump_buffer(std::ofstream& stream, T* start, T* end) { int main() { // const std::vector meshes = {"../../../../mesh/infinite_rectangular_2x8.x"}; - const std::vector meshes = {"../../../../mesh/rectangular_5x10.x"}; + const std::vector meshes = {"../../../../mesh/infinite_rectangular_5x20.x"}; const std::vector backends = get_available_backends(); @@ -77,7 +77,7 @@ int main() { const f32 t_final = 30.0f; const f32 u_inf = 1.0f; // freestream velocity const f32 amplitude = 0.1f; // amplitude of the wing motion - const f32 k = 1.5f; // reduced frequency + const f32 k = 0.5f; // reduced frequency const f32 omega = k * u_inf / (2*b);