Skip to content

Commit

Permalink
better theodorsen plotting
Browse files Browse the repository at this point in the history
  • Loading branch information
@samayala22
samayala22 committed May 21, 2024
1 parent faf8649 commit b784646
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 85 deletions.
102 changes: 49 additions & 53 deletions data/theodorsen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -38,56 +56,45 @@ 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
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)
)

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)
Expand All @@ -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)")

Expand Down
68 changes: 36 additions & 32 deletions tests/test_uvlm_theodorsen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
Expand Down Expand Up @@ -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;
}
}

Expand Down

0 comments on commit b784646

Please sign in to comment.