Skip to content

Commit

Permalink
fix steady force computation
Browse files Browse the repository at this point in the history 
fixes symmetry but scale is still off
  • Loading branch information
@samayala22
samayala22 committed Apr 10, 2024
1 parent 6fbb627 commit 1bbd65f
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 18 deletions.
6 changes: 4 additions & 2 deletions data/theodorsen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,9 @@ 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
b = 0.5 # half chord
a = 10 # full span
a = ar / (2*b) # full span

amplitudes = [0.1, 0.1, 0.1]
reduced_frequencies = [0.5, 0.75, 1.5]
Expand Down Expand Up @@ -76,9 +77,10 @@ def cl_theodorsen(t: float):
uvlm_cl.append(float(cl))

axs["time"].plot(uvlm_t, uvlm_cl, label="UVLM (k=0.5)", linestyle="--")

axs["time"].plot(uvlm_t, uvlm_z, label="UVLM z (k=0.5)", linestyle="--")

axs["heave"].plot(uvlm_z[len(uvlm_cl)//2:], uvlm_cl[len(uvlm_cl)//2:], label="UVLM (k=0.5)", linestyle="--")

axs["time"].set_xlabel('t')
axs["time"].set_ylabel('CL')
axs["time"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray')
Expand Down
53 changes: 53 additions & 0 deletions mesh/rectangular_5x10.x
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
1
6 11 1
0.0000000000000000 0.19999999999999996 0.39999999999999991 0.59999999999999998
0.80000000000000004 1.0000000000000000 0.0000000000000000 0.19999999999999998
0.39999999999999991 0.59999999999999998 0.80000000000000016 1.0000000000000000
0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998
0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000004
0.39999999999999991 0.59999999999999998 0.79999999999999993 1.0000000000000000
0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998
0.80000000000000027 1.0000000000000000 0.0000000000000000 0.20000000000000007
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.80000000000000004 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.20000000000000018 0.39999999999999991 0.60000000000000009
0.79999999999999982 1.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.40000000000000002 0.39999999999999986
0.39999999999999947 0.39999999999999913 0.39999999999999897 0.39999999999999858
0.80000000000000004 0.79999999999999971 0.79999999999999893 0.79999999999999827
0.79999999999999793 0.79999999999999716 1.2000000000000002 1.2000000000000004
1.1999999999999997 1.1999999999999997 1.1999999999999991 1.1999999999999993
1.6000000000000001 1.6000000000000003 1.6000000000000005 1.6000000000000010
1.6000000000000012 1.6000000000000014 2.0000000000000000 2.0000000000000000
2.0000000000000000 2.0000000000000000 2.0000000000000000 2.0000000000000000
2.4000000000000004 2.3999999999999999 2.3999999999999995 2.3999999999999995
2.3999999999999981 2.3999999999999986 2.8000000000000003 2.8000000000000003
2.8000000000000012 2.8000000000000007 2.8000000000000007 2.8000000000000007
3.2000000000000002 3.2000000000000002 3.1999999999999997 3.1999999999999997
3.1999999999999993 3.1999999999999993 3.5999999999999996 3.5999999999999992
3.6000000000000005 3.6000000000000005 3.6000000000000010 3.6000000000000014
4.0000000000000000 4.0000000000000000 4.0000000000000000 4.0000000000000000
4.0000000000000000 4.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
7 changes: 4 additions & 3 deletions tests/test_uvlm_theodorsen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,21 +64,20 @@ void dump_buffer(std::ofstream& stream, T* start, T* end) {

int main() {
// const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x200.x"};
const std::vector<std::string> meshes = {"../../../../mesh/rectangular_5x10.x"};

const std::vector<std::string> backends = get_available_backends();

auto solvers = tiny::make_combination(meshes, backends);


// Geometry
const f32 b = 0.5f; // half chord

// Define simulation length
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 = 0.5f; // reduced frequency
const f32 k = 1.5f; // reduced frequency

const f32 omega = k * u_inf / (2*b);

Expand Down Expand Up @@ -165,6 +164,8 @@ int main() {
std::cout << "position: " << base_vertex << "\n";

const FlowData flow{linalg::alias::float3{-base_velocity[0], -base_velocity[1], -base_velocity[2]}, 1.0f};
std::cout << "lift axis: " << flow.lift_axis << "\n";

backend->compute_rhs(flow);
backend->add_wake_influence(flow);
backend->lu_solve();
Expand Down
64 changes: 51 additions & 13 deletions vlm/backends/cpu/src/vlm_backend_cpu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ void BackendCPU::compute_lhs(const FlowData& flow) {

init.precede(wing_pass, cond);
wing_pass.precede(sync);
cond.precede(wake_pass, sync);
cond.precede(wake_pass, sync); // 0 and 1
wake_pass.precede(back);
back.precede(cond);

Expand Down Expand Up @@ -125,18 +125,49 @@ void BackendCPU::add_wake_influence(const FlowData& flow) {
};

// loop over wake rows
for (u64 i = 0; i < mesh.current_nw; i++) {
const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns;
std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out
// Actually fill the wake buffer
// parallel for
for (u64 j = 0; j < mesh.ns; j++) { // loop over columns
const u64 lidx = wake_row_start + j; // TODO: replace this ASAP
ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
}
// for (u64 i = 0; i < mesh.current_nw; i++) {
// const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns;
// std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out
// // Actually fill the wake buffer
// // parallel for
// for (u64 j = 0; j < mesh.ns; j++) { // loop over columns
// const u64 lidx = wake_row_start + j; // TODO: replace this ASAP
// ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
// }

// cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
// }

u64 idx = 0;
u64 wake_row_start = (m.nc + m.nw - 1) * m.ns;

auto init = taskflow.placeholder();
auto cond = taskflow.emplace([&]{
return idx < mesh.current_nw ? 0 : 1; // 0 means continue, 1 means break
});
auto zero_buffer = taskflow.emplace([&]{
std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out
});
auto wake_influence = taskflow.for_each_index((u64)0, m.ns, [&] (u64 j) {
const u64 lidx = wake_row_start + j;
ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
});
auto back = taskflow.emplace([&]{
cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
}

idx++;
wake_row_start -= m.ns;
return 0; // 0 means continue
});
auto sync = taskflow.placeholder();

init.precede(cond);
cond.precede(zero_buffer, sync);
zero_buffer.precede(wake_influence);
wake_influence.precede(back);
back.precede(cond);

Executor::get().run(taskflow).wait();
}

void BackendCPU::shed_gamma() {
Expand Down Expand Up @@ -216,9 +247,16 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, co
// Steady part:
const linalg::alias::float3 v0 = mesh.get_v0(li);
const linalg::alias::float3 v1 = mesh.get_v1(li);
const linalg::alias::float3 v2 = mesh.get_v2(li);
const linalg::alias::float3 v3 = mesh.get_v3(li);

linalg::alias::float3 force_steady = {0.0f, 0.0f, 0.0f};
force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v1 - v0);
force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v2 - v1);
force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v3 - v2);
force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v0 - v3);

// Leading edge vector pointing outward from wing root
const linalg::alias::float3 dl = v1 - v0;
const linalg::alias::float3 force_steady = flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, dl);
cl += linalg::dot(force_steady, flow.lift_axis);

// Unsteady part (Simpson method)
Expand Down

0 comments on commit 1bbd65f

Please sign in to comment.