From 832b185807005c0cb27101c31f2eaa2e7630658d Mon Sep 17 00:00:00 2001 From: Samuel Ayala Date: Mon, 7 Oct 2024 23:56:07 -0400 Subject: [PATCH] monolithic & iterative dimensionless 2dof solver --- data/2dof.py | 409 +++++++++++++++++++++------------- data/theodorsen_analytical.py | 150 +++++++++++++ 2 files changed, 407 insertions(+), 152 deletions(-) create mode 100644 data/theodorsen_analytical.py diff --git a/data/2dof.py b/data/2dof.py index ef1488a..781b75e 100644 --- a/data/2dof.py +++ b/data/2dof.py @@ -1,7 +1,7 @@ import numpy as np import matplotlib.pyplot as plt -import scipy.integrate as spi -import scipy.special as scsp +from scipy.fft import fft, fftfreq +from scipy.integrate import solve_ivp from tqdm import tqdm EPS_sqrt_f = np.sqrt(1.19209e-07) @@ -42,20 +42,6 @@ def rk2_step(f: callable, x, t, dt, i): k2 = f(t_mid, x_mid) x[i+1] = x[i] + 0.5 * (k1 + k2) * dt -def theo_fun(k): - H1 = scsp.hankel2(1, k) - H0 = scsp.hankel2(0, k) - C = H1 / (H1 + 1.0j * H0) - return C - -# Some info in Katz Plotkin p414 (eq 13.73a) -# Jone's approximation of Wagner function -b0 = 1 -b1 = -0.165 -b2 = -0.335 -beta_1 = 0.0455 -beta_2 = 0.3 - # Dimensional variables class DVars: def __init__(self, rho, u_inf, b, p_axis, x_a, r_a, m, k_h, k_a): @@ -77,151 +63,270 @@ def __init__(self, rho, u_inf, b, p_axis, x_a, r_a, m, k_h, k_a): class NDVars: @classmethod def from_dvars(self, d_vars): - self.p_axis = d_vars.p_axis / d_vars.b # dimensionless pitch axis - self.omega_a = np.sqrt(d_vars.k_a / (d_vars.m * d_vars.r_a)) # dimensionless torsional stiffness - self.omega_h = np.sqrt(d_vars.k_h / d_vars.m) # dimensionless linear stiffness self.U = d_vars.u_inf / (self.omega_a * d_vars.b) # reduced velocity self.mu = d_vars.m / (np.pi * d_vars.rho * d_vars.b**2) # mass ratio self.r_a = d_vars.r_a / d_vars.b # dimensionless radius of gyration self.x_a = d_vars.x_a / d_vars.b # dimensionless mass center - self.eta_omega = self.omega_h / self.omega_a # dimensionless ratio of torsional stiffness to linear stiffness - def __init__(self, a_h, omega, x_a, mu, r_a, U): + def __init__(self, a_h, omega, zeta_a, zeta_h,x_a, mu, r_a, U): self.a_h = a_h # dimensionless pitch axis self.omega = omega # dimensionless ratio of torsional stiffness to linear stiffness + self.zeta_a = zeta_a # dimensionless pitch axis + self.zeta_h = zeta_h # dimensionless pitch axis self.x_a = x_a # dimensionless mass center self.mu = mu # mass ratio self.r_a = r_a # dimensionless radius of gyration self.U = U # reduced velocity -d_vars = DVars( - rho = 1, # fluid density - u_inf = 1, # freestream - b = 0.5, # half chord - p_axis = -0.25, # pitch/elastic axis (measured from center of wing) - x_a = 0.125, # mass center (from pitch axis to center of mass) - r_a = 0.25, # radius of gyration around elastic axis - m = 1.0, # mass - k_h = 800.0, # linear stiffness - k_a = 150.0 # torsional stiffness -) - -# Zhao 2004 params -nd_vars = NDVars( - a_h = -0.5, - omega = 0.2, - x_a = 0.25, - mu = 100.0, - r_a = 0.5 -) - -# Theodorsen numerical simulation param -dt = 0.01 -t_final = 30 # 30s -vec_t = np.arange(0, t_final, dt) -vec_t_nd = vec_t * -# Aeroelastic model -# m = 2.0 # mass -# S_a = 0.2 # first moment of inertia -# I_a = 4.0 # second moment of inertia -# K_h = 800.0 # linear stiffness -# K_a = 150.0 # torsional stiffness - -M = np.array([ - [1/(nd_vars.r_a**2), nd_vars.x_a/(nd_vars.r_a**2)], - [nd_vars.x_a, 1.0] -]) -print(M) -C = np.zeros((2,2)) -K = np.array([[nd_vars.eta_omega, 0],[0, 1]]) - -u = np.zeros((2, len(vec_t))) -v = np.zeros((2, len(vec_t))) -a = np.zeros((2, len(vec_t))) -F = np.zeros((2, len(vec_t))) - -# augmented aero states -vec_x1 = np.zeros(len(vec_t)) -vec_x2 = np.zeros(len(vec_t)) - -# Initial condition -u[:, 0] = np.array([0.0, np.radians(3)]) -v[:, 0] = np.array([0, 0]) -a[:, 0] = np.linalg.solve(M, F[:, 0] - C @ v[:,0] - K @ u[:,0]) - -fig, axs = plt.subplot_mosaic( - [["h"], ["alpha"]], # Disposition des graphiques - constrained_layout=True, # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas - figsize=(11, 8), # Ajuster la taille de la figure (x,y) -) - -def w(s: float): - idx = int(s / dt) - return d_vars.u_inf * u[1, idx] + v[0,idx] + d_vars.b * (0.5 - nd_vars.p_axis) * v[1, idx] - -def dx1ds(s: float, x1: float): return (d_vars.u_inf / d_vars.b) * (b1 * beta_1 * w(s) - beta_1 * x1) -def dx2ds(s: float, x2: float): return (d_vars.u_inf / d_vars.b) * (b2 * beta_2 * w(s) - beta_2 * x2) - -def aero(i): - t = vec_t[i] - - rk2_step(dx1ds, vec_x1, vec_t[i-1], dt, i-1) - rk2_step(dx2ds, vec_x2, vec_t[i-1], dt, i-1) - x1 = vec_x1[i] - x2 = vec_x2[i] - - L_m = d_vars.rho * d_vars.b * d_vars.b * np.pi * (d_vars.u_inf * v[1,i] + a[0,i] - d_vars.b * nd_vars.p_axis * a[1,i]) - L_c = -2 * np.pi * d_vars.rho * d_vars.u_inf * d_vars.b * (-(b0 + b1 + b2) * w(t) + x1 + x2) - L = L_m + L_c - M = (nd_vars.p_axis + 0.5) * d_vars.b * np.cos(u[1,i]) * L - F[0,i] = -L - F[1,i] = M - # C_l = L / (0.5 * rho * u_inf * u_inf * c) - -# def central_difference_step(M, C, K, u_prev, u, F_i, dt): -# M_inv = np.linalg.inv(M) -# a = M_inv @ (F_i - C @ ((u - u_prev) / dt) - K @ u) -# u_next = 2 * u - u_prev + dt**2 * a -# return u_next - -# # Initialize u_prev -# u_prev = u[:,0] - dt * v[:,0] + 0.5 * dt**2 * a[:,0] -# for i in tqdm(range(1, len(vec_t))): -# F_i = F[:,i-1] -# u[:,i] = central_difference_step(M, C, K, u_prev, u[:,i-1], F_i, dt) -# v[:,i] = (u[:,i] - u_prev) / (2 * dt) -# a[:,i] = (u[:,i] - 2 * u[:,i-1] + u_prev) / dt**2 -# u_prev = u[:,i-1] -# aero(i) - -# Newmark V2 -for i in tqdm(range(len(vec_t)-1)): - t = vec_t[i] - F[:,i+1] = F[:,i] - du = np.zeros(2) - du_k = np.zeros(2) + 1 - iteration = 0 - while (np.linalg.norm(du_k - du) / len(du) > 1e-7): - du_k = du[:] - du, dv, da = newmark_beta_step(M, C, K, u[:,i], v[:,i], a[:,i], F[:,i], F[:,i+1], dt) - u[:,i+1] = u[:,i] + du - v[:,i+1] = v[:,i] + dv - a[:,i+1] = a[:,i] + da - aero(i+1) - iteration += 1 - # print("iters: ", iteration) - -axs["h"].plot(vec_t, u[0, :]) -axs["alpha"].plot(vec_t, u[1, :]) - -axs["h"].set_xlabel('t') -axs["h"].set_ylabel('h') -axs["h"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') - -axs["alpha"].set_xlabel('t') -axs["alpha"].set_ylabel('alpha') -axs["alpha"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') - -# plt.suptitle("Verification of UVLM with Theodorsen harmonic heave motion") - -plt.show() +def monolithic_aeroelastic_system(ndv: NDVars): + psi1 = 0.165 + psi2 = 0.335 + eps1 = 0.0455 + eps2 = 0.3 + + c0 = 1 + 1 / ndv.mu + c1 = ndv.x_a - ndv.a_h / ndv.mu + c2 = (2 / ndv.mu) * (1 - psi1 - psi2) + c3 = (1 / ndv.mu) * (1 + (1 - 2 * ndv.a_h) * (1 - psi1 - psi2)) + c4 = (2 / ndv.mu) * (eps1 * psi1 + eps2 * psi2) + c5 = (2 / ndv.mu) * (1 - psi1 - psi2 + (0.5 - ndv.a_h) * (eps1 * psi1 + eps2 * psi2)) + c6 = (2 / ndv.mu) * eps1 * psi1 * (1 - eps1 * (0.5 - ndv.a_h)) + c7 = (2 / ndv.mu) * eps2 * psi2 * (1 - eps2 * (0.5 - ndv.a_h)) + c8 = -(2 / ndv.mu) * eps1**2 * psi1 + c9 = -(2 / ndv.mu) * eps2**2 * psi2 + + r_alpha_sq = ndv.r_a ** 2 + d0 = (ndv.x_a / r_alpha_sq) - (ndv.a_h / (ndv.mu * r_alpha_sq)) + d1 = 1 + (1 + 8 * ndv.a_h**2) / (8 * ndv.mu * r_alpha_sq) + d2 = -((1 + 2 * ndv.a_h) / (ndv.mu * r_alpha_sq)) * (1 - psi1 - psi2) + d3 = ((1 - 2 * ndv.a_h) / (2 * ndv.mu * r_alpha_sq)) - \ + ((1 + 2 * ndv.a_h) * (1 - 2 * ndv.a_h) * (1 - psi1 - psi2)) / (2 * ndv.mu * r_alpha_sq) + d4 = -((1 + 2 * ndv.a_h) / (ndv.mu * r_alpha_sq)) * (eps1 * psi1 + eps2 * psi2) + d5 = -((1 + 2 * ndv.a_h) / (ndv.mu * r_alpha_sq)) * (1 - psi1 - psi2) - \ + ((1 + 2 * ndv.a_h) * (1 - 2 * ndv.a_h) * (psi1 * eps1 - psi2 * eps2)) / (2 * ndv.mu * r_alpha_sq) + d6 = -((1 + 2 * ndv.a_h) * psi1 * eps1) / (ndv.mu * r_alpha_sq) * (1 - eps1 * (0.5 - ndv.a_h)) + d7 = -((1 + 2 * ndv.a_h) * psi2 * eps2) / (ndv.mu * r_alpha_sq) * (1 - eps2 * (0.5 - ndv.a_h)) + d8 = ((1 + 2 * ndv.a_h) * psi1 * eps1**2) / (ndv.mu * r_alpha_sq) + d9 = ((1 + 2 * ndv.a_h) * psi2 * eps2**2) / (ndv.mu * r_alpha_sq) + + j = 1 / (c0 * d1 - c1 * d0) + a21 = j * (-d5 * c0 + c5 * d0) + a22 = j * (-d3 * c0 + c3 * d0) + a23 = j * (-d4 * c0 + c4 * d0) + a24 = j * (-d2 * c0 + c2 * d0) + a25 = j * (-d6 * c0 + c6 * d0) + a26 = j * (-d7 * c0 + c7 * d0) + a27 = j * (-d8 * c0 + c8 * d0) + a28 = j * (-d9 * c0 + c9 * d0) + a41 = j * (d5 * c1 - c5 * d1) + a42 = j * (d3 * c1 - c3 * d1) + a43 = j * (d4 * c1 - c4 * d1) + a44 = j * (d2 * c1 - c2 * d1) + a45 = j * (d6 * c1 - c6 * d1) + a46 = j * (d7 * c1 - c7 * d1) + a47 = j * (d8 * c1 - c8 * d1) + a48 = j * (d9 * c1 - c9 * d1) + + A = np.zeros((8, 8)) + A[0, 1] = 1 + A[1, 0] = a21 - j*c0*(1/ndv.U)**2 + A[1, 1] = a22 + 2*j*c0*ndv.zeta_a*(1/ndv.U) + A[1, 2] = a23 + j*d0*(ndv.omega/ndv.U)**2 + A[1, 3] = a24 + 2*j*d0*ndv.zeta_h*(ndv.omega/ndv.U) + A[1, 4] = a25 + A[1, 5] = a26 + A[1, 6] = a27 + A[1, 7] = a28 + A[2, 3] = 1 + A[3, 0] = a41 + j*c1*(1/ndv.U)**2 + A[3, 1] = a42 + 2*j*c1*ndv.zeta_a*(1/ndv.U) + A[3, 2] = a43 - j*d1*(ndv.omega/ndv.U)**2 + A[3, 3] = a44 - 2*j*d1*ndv.zeta_h*(ndv.omega/ndv.U) + A[3, 4] = a45 + A[3, 5] = a46 + A[3, 6] = a47 + A[3, 7] = a48 + A[4, 0] = 1 + A[4, 4] = -eps1 + A[5, 0] = 1 + A[5, 5] = -eps2 + A[6, 2] = 1 + A[6, 6] = -eps1 + A[7, 3] = 1 + A[7, 7] = -eps2 + + return A + +# d_vars = DVars( +# rho = 1, # fluid density +# u_inf = 1, # freestream +# b = 0.5, # half chord +# p_axis = -0.25, # pitch/elastic axis (measured from center of wing) +# x_a = 0.125, # mass center (from pitch axis to center of mass) +# r_a = 0.25, # radius of gyration around elastic axis +# m = 1.0, # mass +# k_h = 800.0, # linear stiffness +# k_a = 150.0 # torsional stiffness +# ) + +if __name__ == "__main__": + psi1 = 0.165 + psi2 = 0.335 + eps1 = 0.0455 + eps2 = 0.3 + + # vec_U = np.linspace(0.1, 7, 50) + vec_U = [3.0] + freqs = np.zeros((2, len(vec_U))) + + # Dimensionless parameters + dt_nd = 0.01 + t_final_nd = 300 + vec_t_nd = np.arange(0, t_final_nd, dt_nd) + n = len(vec_t_nd) + + for u_idx, U_vel in enumerate(vec_U): + # Zhao 2004 params + ndv = NDVars( + a_h = -0.5, + omega = 0.2, + zeta_a = 0.0, + zeta_h = 0.0, + x_a = 0.25, + mu = 100.0, + r_a = 0.5, + U = U_vel + ) + + # Aeroelastic equations of motion mass, damping and stiffness matrices + M = np.array([ + [1.0, ndv.x_a], + [ndv.x_a / (ndv.r_a**2), 1.0] + ]) + C = np.array([ + [2.0*ndv.zeta_h*(ndv.omega/ndv.U), 0], + [0, 2.0*ndv.zeta_a/ndv.U] + ]) + K = np.array([ + [(ndv.omega/ndv.U)**2, 0], + [0, 1/(ndv.U**2)] + ]) + + # Time history of dofs + u = np.zeros((2, len(vec_t_nd))) + v = np.zeros((2, len(vec_t_nd))) + a = np.zeros((2, len(vec_t_nd))) + F = np.zeros((2, len(vec_t_nd))) + + # Augmented aero states + vec_w1 = np.zeros(len(vec_t_nd)) + vec_w2 = np.zeros(len(vec_t_nd)) + + # Initial condition + u[:, 0] = np.array([0.0, np.radians(3)]) + v[:, 0] = np.array([0, 0]) + a[:, 0] = np.linalg.solve(M, F[:, 0] - C @ v[:,0] - K @ u[:,0]) + + def w(s: float): + idx = int(s / dt_nd) + return v[1, idx] + a[0, idx] + (0.5 - ndv.a_h) * a[1, idx] + + def dw1ds(s: float, w1: float): return w(s) - eps1 * w1 + def dw2ds(s: float, w2: float): return w(s) - eps2 * w2 + + def aero(i): + t = vec_t_nd[i] + + rk2_step(dw1ds, vec_w1, vec_t_nd[i-1], dt_nd, i-1) + rk2_step(dw2ds, vec_w2, vec_t_nd[i-1], dt_nd, i-1) + w1 = vec_w1[i] + w2 = vec_w2[i] + + duhamel = u[1, i] - u[1, 0] + v[0, i] - v[0, 0] + (0.5 - ndv.a_h)*(v[1, i] - v[1, 0]) - psi1*w1 - psi2*w2 + wagner_init = (u[1, 0] + u[0, 0] + (0.5 - ndv.a_h)*v[1,0])*(1 - psi1*np.exp(-eps1*t) - psi2*np.exp(-eps2*t)) + cl = np.pi*(a[0, i] - ndv.a_h * a[1, i] + v[1, i]) + 2*np.pi*(wagner_init + duhamel) + cm = np.pi*(0.5 + ndv.a_h)*(wagner_init + duhamel) + 0.5*np.pi*ndv.a_h*(a[0, i] - ndv.a_h*a[1, i]) - 0.5*np.pi*(0.5 - ndv.a_h)*v[1, i] - (np.pi/16) * a[1, i] + F[0,i] = - cl / (np.pi*ndv.mu) + F[1,i] = (2*cm) / (np.pi*ndv.mu*ndv.r_a**2) + + # def central_difference_step(M, C, K, u_prev, u, F_i, dt): + # M_inv = np.linalg.inv(M) + # a = M_inv @ (F_i - C @ ((u - u_prev) / dt) - K @ u) + # u_next = 2 * u - u_prev + dt**2 * a + # return u_next + + # # Initialize u_prev + # u_prev = u[:,0] - dt * v[:,0] + 0.5 * dt**2 * a[:,0] + # for i in tqdm(range(1, len(vec_t))): + # F_i = F[:,i-1] + # u[:,i] = central_difference_step(M, C, K, u_prev, u[:,i-1], F_i, dt) + # v[:,i] = (u[:,i] - u_prev) / (2 * dt) + # a[:,i] = (u[:,i] - 2 * u[:,i-1] + u_prev) / dt**2 + # u_prev = u[:,i-1] + # aero(i) + + # Newmark V2 + for i in tqdm(range(n-1)): + t = vec_t_nd[i] + F[:,i+1] = F[:,i] + du = np.zeros(2) + du_k = np.zeros(2) + 1 + iteration = 0 + while (np.linalg.norm(du_k - du) / len(du) > 1e-7): + du_k = du[:] + du, dv, da = newmark_beta_step(M, C, K, u[:,i], v[:,i], a[:,i], F[:,i], F[:,i+1], dt_nd) + u[:,i+1] = u[:,i] + du + v[:,i+1] = v[:,i] + dv + a[:,i+1] = a[:,i] + da + aero(i+1) + iteration += 1 + # print("iters: ", iteration) + + A = monolithic_aeroelastic_system(ndv) + y0 = np.array([np.radians(3), 0, 0, 0, 0, 0, 0, 0]) + monolithic_sol = solve_ivp(lambda t, y: A @ y, (0, t_final_nd), y0, t_eval=vec_t_nd) + + fig, axs = plt.subplot_mosaic( + [["h"], ["alpha"]], # Disposition des graphiques + constrained_layout=True, # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas + figsize=(11, 8), # Ajuster la taille de la figure (x,y) + ) + axs["h"].plot(vec_t_nd, u[0, :], label="Iterative") + axs["alpha"].plot(vec_t_nd, u[1, :], label="Iterative") + axs["h"].plot(vec_t_nd, monolithic_sol.y[2, :], label="Monolithic") + axs["alpha"].plot(vec_t_nd, monolithic_sol.y[0, :], label="Monolithic") + + axs["h"].legend() + axs["alpha"].legend() + axs["h"].set_xlabel('t') + axs["h"].set_ylabel('h') + axs["h"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') + + axs["alpha"].set_xlabel('t') + axs["alpha"].set_ylabel('alpha') + axs["alpha"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') + plt.show() + + # X = fft(u, axis=1) + # freq = fftfreq(n, d=dt_nd) + # idx = np.argmax(np.abs(X), axis=1) + # dominant_freqs = np.abs(freq[idx]) + # freqs[:, u_idx] = dominant_freqs + + # fig, axs = plt.subplot_mosaic( + # [["freq"]], # Disposition des graphiques + # constrained_layout=True, # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas + # figsize=(11, 8), # Ajuster la taille de la figure (x,y) + # ) + # axs["freq"].plot(vec_U, freqs[0, :], "o") + # axs["freq"].plot(vec_U, freqs[1, :], "o") + # plt.show() + + # fig, axs = plt.subplot_mosaic( + # [["damping"]], # Disposition des graphiques + # constrained_layout=True, # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas + # figsize=(11, 8), # Ajuster la taille de la figure (x,y) + # ) + # axs["damping"].plot(vec_U, freqs[0, :], "o") + # axs["damping"].plot(vec_U, freqs[1, :], "o") + # plt.show() diff --git a/data/theodorsen_analytical.py b/data/theodorsen_analytical.py new file mode 100644 index 0000000..bd09653 --- /dev/null +++ b/data/theodorsen_analytical.py @@ -0,0 +1,150 @@ +import numpy as np +import matplotlib.pyplot as plt +import scipy.integrate as spi + +EPS_sqrt_f = np.sqrt(1.19209e-07) + +# Two point central difference +def derivative(f, x): + return (f(x + EPS_sqrt_f) - f(x - EPS_sqrt_f)) / (2 * EPS_sqrt_f) + +def derivative2(f, x): + return (f(x + EPS_sqrt_f) - 2 * f(x) + f(x - EPS_sqrt_f)) / (EPS_sqrt_f ** 2) + +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 cl_theodorsen1(t: np.ndarray, alpha:callable, h:callable, a_h: float, b: float): + """ + t: time + alpha: angle of attack + h: heave + a_h: dimensionless distance from center chord to pitch axis + returns cl vector + """ + # Jone's approximation of Wagner function + b0 = 1 + b1 = -0.165 + b2 = -0.335 + beta_1 = 0.0455 + beta_2 = 0.3 + + def w(s: float): return u_inf * alpha(s) + derivative(h, s) + b * (0.5 - a_h) * derivative(alpha, s) + def dx1ds(s: float, x1: float): return (u_inf / b) * (b1 * beta_1 * w(s) - beta_1 * x1) + def dx2ds(s: float, x2: float): return (u_inf / b) * (b2 * beta_2 * w(s) - beta_2 * x2) + x1_solution = spi.solve_ivp(dx1ds, [0, t[-1]], [0], t_eval=t) + x2_solution = spi.solve_ivp(dx2ds, [0, t[-1]], [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]) + L_m = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * a_h * 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) + +# 4 augmented variables +def cl_theodorsen2(t: np.ndarray, alpha:callable, h:callable, a_h: float): + """ + t: dimensionlesstime + alpha: angle of attack + h: dimensionlessheave + a_h: dimensionless distance from center chord to pitch axis + returns cl vector + """ + psi1 = 0.165 + psi2 = 0.335 + eps1 = 0.0455 + eps2 = 0.3 + def dw1ds(s: float, w1: float): return alpha(s) - eps1 * w1 + def dw2ds(s: float, w2: float): return alpha(s) - eps2 * w2 + def dw3ds(s: float, w3: float): return h(s) - eps1 * w3 + def dw4ds(s: float, w4: float): return h(s) - eps2 * w4 + w1_solution = spi.solve_ivp(dw1ds, [0, t[-1]], [0], t_eval=t) + w2_solution = spi.solve_ivp(dw2ds, [0, t[-1]], [0], t_eval=t) + w3_solution = spi.solve_ivp(dw3ds, [0, t[-1]], [0], t_eval=t) + w4_solution = spi.solve_ivp(dw4ds, [0, t[-1]], [0], t_eval=t) + def w1(s: float): return np.interp(s, w1_solution.t, w1_solution.y[0]) + def w2(s: float): return np.interp(s, w2_solution.t, w2_solution.y[0]) + def w3(s: float): return np.interp(s, w3_solution.t, w3_solution.y[0]) + def w4(s: float): return np.interp(s, w4_solution.t, w4_solution.y[0]) + part0 = (alpha(t) - alpha(0) + derivative(h, t) - derivative(h, 0) + (0.5 - a_h)*(derivative(alpha, t) - derivative(alpha, 0))) + part11 = -psi1*(alpha(t) - np.exp(-eps1*t)*alpha(0) - eps1*w1(t)) + part12 = -psi1*(derivative(h, t) - np.exp(-eps1*t)*derivative(h, 0) - eps1*h(t) + eps1 * np.exp(-eps1*t)*h(0) +eps1**2*w3(t)) + part13 = -psi1*(0.5-a_h)*(derivative(alpha, t) - np.exp(-eps1*t)*derivative(alpha, 0) - eps1*alpha(t) + eps1 * np.exp(-eps1*t)*alpha(0) + eps1**2*w1(t)) + part21 = -psi2*(alpha(t) - np.exp(-eps2*t)*alpha(0) - eps2*w2(t)) + part22 = -psi2*(derivative(h, t) - np.exp(-eps2*t)*derivative(h, 0) - eps2*h(t) + eps2 * np.exp(-eps2*t)*h(0) +eps2**2*w4(t)) + part23 = -psi2*(0.5-a_h)*(derivative(alpha, t) - np.exp(-eps2*t)*derivative(alpha, 0) - eps2*alpha(t) + eps2 * np.exp(-eps2*t)*alpha(0) + eps2**2*w2(t)) + cl = np.pi*(derivative2(h, t) - a_h * derivative2(alpha, t) + derivative(alpha, t)) + 2*np.pi*(alpha(0) + derivative(h, 0) + (0.5 - a_h)*derivative(alpha, 0))*(1-psi1*np.exp(-eps1*t)-psi2*np.exp(-eps2*t))+2*np.pi*(part0+part11+part12+part13+part21+part22+part23) + return cl + +# 2 augmented variables +def cl_theodorsen3(t: np.ndarray, alpha:callable, h:callable, a_h: float): + """ + t: dimensionlesstime + alpha: angle of attack + h: dimensionlessheave + a_h: dimensionless distance from center chord to pitch axis + returns cl vector + """ + psi1 = 0.165 + psi2 = 0.335 + eps1 = 0.0455 + eps2 = 0.3 + def w(s: float): return derivative(alpha, s) + derivative2(h, s) + (0.5 - a_h) * derivative2(alpha, s) + def dw1ds(s: float, w1: float): return w(s) - eps1 * w1 + def dw2ds(s: float, w2: float): return w(s) - eps2 * w2 + w1_solution = spi.solve_ivp(dw1ds, [0, t[-1]], [0], t_eval=t) + w2_solution = spi.solve_ivp(dw2ds, [0, t[-1]], [0], t_eval=t) + def w1(s: float): return np.interp(s, w1_solution.t, w1_solution.y[0]) + def w2(s: float): return np.interp(s, w2_solution.t, w2_solution.y[0]) + part0 = (alpha(t) + derivative(h, t) + (0.5 - a_h)*derivative(alpha, t)) + part1 = -psi1*w1(t) + part2 = -psi2*w2(t) + cl = np.pi*(derivative2(h, t) - a_h * derivative2(alpha, t) + derivative(alpha, t)) + 2*np.pi*(alpha(0) + derivative(h, 0) + (0.5 - a_h)*derivative(alpha, 0))*(-psi1*np.exp(-eps1*t)-psi2*np.exp(-eps2*t))+2*np.pi*(part0+part1+part2) + return cl + +if __name__ == "__main__": + # UVLM parameters + rho = 1 # fluid density + u_inf = 1 # freestream + ar = 10000 # aspect ratio + b = 0.5 # half chord + c = 2*b # chord + a = ar / c # full span + a_h = -0.5 # quarter chord + + nb_pts = 1000 + k = 0.5 # reduced frequency + t_final = 60 # dimensional time + omega = k * 2.0 * u_inf / c # frequency + vec_t = np.linspace(0, t_final, nb_pts) + + # 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 -0.1 * np.sin(omega * t) + + # pure pitching + def pitch(t): return np.radians(np.sin(omega * t)) + def heave(t): return 0 + + fig, axs = plt.subplot_mosaic( + [["time"]], # Disposition des graphiques + constrained_layout=True, # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas + figsize=(11, 6), # Ajuster la taille de la figure (x,y) + ) + + axs["time"].plot(vec_t, cl_theodorsen1(vec_t, pitch, heave, a_h, b), label=f"Theodorsen (k={k})") + + vec_t_nd = u_inf * vec_t / b # non-dimensional time + axs["time"].plot(vec_t, cl_theodorsen3(vec_t_nd, lambda t: pitch(t*b / u_inf), lambda t: heave(t*b / u_inf) / b, a_h), label=f"Theodorsen 2 (k={k})") + + axs["time"].set_xlabel('t') + axs["time"].set_ylabel('CL') + axs["time"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') + axs["time"].legend() + + plt.suptitle("Theodorsen solution") + + plt.show()