From ac067532b3c58729a54fd502d2757b0ac24f8c06 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Corentin=20Thom=C3=A9e?= <cthomee@uliege.be>
Date: Mon, 20 Jan 2025 16:19:42 +0100
Subject: [PATCH] CutTE
---
hspm/api/core.py | 11 +++++++
hspm/src/geometry.cpp | 70 +++++++++++++++++++++++++++++++++++++++++--
hspm/src/hspm.h | 4 +++
hspm/src/solver.cpp | 8 ++---
4 files changed, 87 insertions(+), 6 deletions(-)
diff --git a/hspm/api/core.py b/hspm/api/core.py
index 7905e39..8787ed3 100644
--- a/hspm/api/core.py
+++ b/hspm/api/core.py
@@ -63,6 +63,17 @@ def initHSPM(cfg):
hspm.sLe = cfg['sLe']
hspm.sTe = cfg['sTe']
cfg["N"] = 1 # Placeholder value for the generateNaca4DigitCoordinates function
+
+ if 'cutTE' in cfg:
+ hspm.cutTE = cfg['cutTE']
+ if not "sTe" in cfg:
+ print("HSPM Warning: Trailing edge mesh size (sTe) not defined. Using default value of 0.01")
+ cfg['sTe'] = 0.01
+ hspm.sTe = cfg['sTe']
+
+ if 'stopTE' in cfg:
+ hspm.stopTE = True
+ hspm.stopTELocation = cfg['stopTE']
if cfg['naca']:
hspm.generateNaca4DigitCoordinates(int(cfg['naca'][0]), int(cfg['naca'][1]), int(cfg['naca'][2:4]), int(cfg['N']))
diff --git a/hspm/src/geometry.cpp b/hspm/src/geometry.cpp
index 6c1321e..315a0d0 100644
--- a/hspm/src/geometry.cpp
+++ b/hspm/src/geometry.cpp
@@ -67,8 +67,6 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
}
N = _N;
- x = Eigen::VectorXd(N+1);
- y = Eigen::VectorXd(N+1);
Eigen::VectorXd angleDistribution = Eigen::VectorXd::LinSpaced(N+1, 0, 2*M_PI);
@@ -117,6 +115,74 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
}
}
+ if (cutTE) {
+ size_t newN = N;
+ size_t remPts = 0;
+ for (size_t i=1; i<=N/2; i++) {
+ if (chord - x_distr(i) >= sTe) {
+ std::cout << "Removed " << i-1 << " points from the trailing edge. " << std::endl;
+ newN = N - (i-1)*2;
+ remPts = i-1;
+ break;
+ }
+ }
+
+ Eigen::VectorXd x_distr_new = Eigen::VectorXd(newN+1);
+ x_distr_new(0) = x_distr(0);
+ for (size_t i=1; i<newN; i++) {
+ x_distr_new(i) = x_distr(i+remPts);
+ }
+ x_distr_new(newN) = x_distr(N);
+
+ N = newN;
+ x_distr = x_distr_new;
+ std::cout << "HSPM: The mesh contains " << N << " elements." << std::endl;
+ }
+
+ if (stopTE) {
+ // First, find the mesh size at stopTELOcation
+ double stopSize;
+ size_t N_behind;
+ double actualStopTELocation;
+ double behindSize;
+ size_t N_rem;
+ size_t newN;
+ for (size_t i=1; i<N/2; i++) {
+ if (x_distr(i) <= stopTELocation) {
+ // We found the location
+ actualStopTELocation = x_distr(i);
+ stopSize = x_distr(i) - x_distr(i+1);
+ N_behind = (1 - actualStopTELocation) / stopSize;
+ behindSize = (1 - actualStopTELocation) / N_behind;
+ N_rem = i;
+
+ newN = N+2*(N_behind-N_rem);
+ break;
+ }
+ }
+
+ std::cout << newN << std::endl;
+
+ Eigen::VectorXd x_distr_new = Eigen::VectorXd(newN+1);
+ for (size_t i=0; i<N_behind; i++) {
+ x_distr_new(i) = 1-behindSize*i;
+ }
+
+ for (size_t i=0; i<N-2*N_rem; i++) {
+ x_distr_new(i+N_behind) = x_distr(i+N_rem);
+ }
+ for (size_t i=0; i<N_behind+1; i++) {
+ x_distr_new(newN-i) = 1 - i*behindSize;
+ }
+
+ N = newN;
+ x_distr = x_distr_new;
+ std::cout << "HSPM: The mesh contains " << N << " elements." << std::endl;
+ }
+
+ x = Eigen::VectorXd(N+1);
+ y = Eigen::VectorXd(N+1);
+
for (size_t i=0; i<=N; i++) {
double T = 10 * thickness/100 * chord * (.2969*sqrt(x_distr(i)/chord) - .126*x_distr(i)/chord
- .3537*pow(x_distr(i)/chord, 2) + .2843*pow(x_distr(i)/chord, 3) - .1015*pow(x_distr(i)/chord, 4));
diff --git a/hspm/src/hspm.h b/hspm/src/hspm.h
index 0e37ce4..73beefd 100644
--- a/hspm/src/hspm.h
+++ b/hspm/src/hspm.h
@@ -84,6 +84,10 @@ public:
double sLe, sTe;
double prog; // For exponential mesh
+ bool cutTE;
+ bool stopTE;
+ double stopTELocation;
+
private:
Eigen::PartialPivLU<Eigen::MatrixXd> solver;
};
diff --git a/hspm/src/solver.cpp b/hspm/src/solver.cpp
index 3cadd30..6039d41 100644
--- a/hspm/src/solver.cpp
+++ b/hspm/src/solver.cpp
@@ -133,7 +133,7 @@ void HSPM::computePressureDistribution() {
// Compute the tangential speeds
for (size_t i=0; i<N; i++) {
- V_t(i) = V_inf * cos(AoA - theta(i));
+ /*V_t(i) = V_inf * cos(AoA - theta(i));
for (size_t j=0; j<N; j++) {
V_t(i) += A_t(i, j) * q(j);
V_t(i) += B_t(i, j) * tau;
@@ -143,9 +143,9 @@ void HSPM::computePressureDistribution() {
Cp(i, 0) = x_m(i);
Cp(i, 1) = y_m(i);
Cp(i, 2) = 0;
- Cp(i, 3) = 1 - pow(V_t(i) / V_inf, 2);
+ Cp(i, 3) = 1 - pow(V_t(i) / V_inf, 2);*/
- /*
+
// Off body Cp, gives better graphs !
double x_visc = x_m(i) - dStar(i) * sin(theta(i));
double y_visc = y_m(i) + dStar(i) * cos(theta(i));
@@ -160,7 +160,7 @@ void HSPM::computePressureDistribution() {
Cp(i, 1) = y_m(i);
Cp(i, 2) = 0;
Cp(i, 3) = 1 - pow(V_mag / V_inf, 2);
- */
+
}
cl = 0;
--
GitLab