diff --git a/hspm/src/geometry.cpp b/hspm/src/geometry.cpp
index 2c54f3eb76cf3e26cd6df585cb99e897ab39588f..e8ec92695b9f5b57fff66056220fdaf654c2bd1d 100644
--- a/hspm/src/geometry.cpp
+++ b/hspm/src/geometry.cpp
@@ -23,13 +23,14 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
Eigen::VectorXd x_distr = Eigen::VectorXd::Zero(N+1);
for (size_t i=0; i<=N; i++) {
double angle = angleDistribution(i);
- if (angle < M_PI/2) {
+ /*if (angle < M_PI/2) {
x_distr(i) = chord/2 * (1 - angle * 2 / M_PI + 1);
} else if (angle > 3*M_PI/2) {
x_distr(i) = chord/2 * ((angle - 3*M_PI/2) * 2 / M_PI + 1);
} else {
x_distr(i) = chord/2 * (cos(angle)+1);
- }
+ }*/
+ x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)),2);
}
for (size_t i=0; i<=N; i++) {
diff --git a/hspm/src/hspm.cpp b/hspm/src/hspm.cpp
index f37c052ec0fcc17ce59f188adabb20befd6b24f7..bd7be9c29f4a52dd6655f247296d610ee4fac2c9 100644
--- a/hspm/src/hspm.cpp
+++ b/hspm/src/hspm.cpp
@@ -13,4 +13,20 @@ void HSPM::initHSPM()
this->dStar = Eigen::VectorXd::Zero(N);
V_x = V_inf * cos(AoA);
V_y = V_inf * sin(AoA);
+
+
+ // TODO: Temporary
+ this->L_wake = 3 * chord;
+ this->N_wake = 100;
+ this->x_wake = Eigen::VectorXd::LinSpaced(N_wake+1, chord, L_wake+chord);
+ this->y_wake = Eigen::VectorXd::Zero(N_wake+1);
+ this->x_m_wake = Eigen::VectorXd::Zero(N_wake);
+ this->y_m_wake = Eigen::VectorXd::Zero(N_wake);
+ for (size_t i=0; i<N_wake; i++) {
+ this->x_m_wake(i) = (x_wake(i) + x_wake(i+1)) / 2;
+ this->y_m_wake(i) = (y_wake(i) + y_wake(i+1)) / 2;
+ }
+ this->U_wake = Eigen::VectorXd::Zero(N_wake);
+ this->V_wake = Eigen::VectorXd::Zero(N_wake);
+
}
\ No newline at end of file
diff --git a/hspm/src/hspm.h b/hspm/src/hspm.h
index d65917c71e4d47bf8924ca1bda60bbe21a089198..ee2507a1ee51fc6381957c2f777023e5a19ae852 100644
--- a/hspm/src/hspm.h
+++ b/hspm/src/hspm.h
@@ -63,6 +63,15 @@ public:
double cl;
double cm;
+ // Wake stuff
+ double L_wake;
+ size_t N_wake;
+ Eigen::VectorXd x_wake;
+ Eigen::VectorXd y_wake;
+ Eigen::VectorXd x_m_wake;
+ Eigen::VectorXd y_m_wake;
+ Eigen::VectorXd U_wake;
+ Eigen::VectorXd V_wake;
private:
};
diff --git a/hspm/src/influenceCoeffs.cpp b/hspm/src/influenceCoeffs.cpp
index 9b263e2c7a96b6bcca90bf801b03a5871746d6e6..3ee327f3e7a2bc2e896e1fa8b708060b1045528f 100644
--- a/hspm/src/influenceCoeffs.cpp
+++ b/hspm/src/influenceCoeffs.cpp
@@ -20,7 +20,7 @@ void HSPM::computeConstantInfluenceCoeffs()
theta(i) = atan2(y(i) - y(i+1), x(i) - x(i+1));
- lengths(i) = sqrt(pow(x((i+1)%N) - x(i), 2) + pow(y((i+1)%N) - y(i), 2));
+ lengths(i) = sqrt(pow(x(i+1) - x(i), 2) + pow(y(i+1) - y(i), 2));
for (size_t j=0; j<N; j++) {
beta(i, j) = atan2(y_m(i) - y(j), x_m(i) - x(j)) - atan2(y_m(i) - y(j+1), x_m(i) - x(j+1));
diff --git a/hspm/src/interface.cpp b/hspm/src/interface.cpp
index 1907ef1064219df284b0c46dc2a58bc96a4efb08..42b2473d3d8d5a150eb0a49f3b3b0bf3ef58fafd 100644
--- a/hspm/src/interface.cpp
+++ b/hspm/src/interface.cpp
@@ -12,7 +12,8 @@ void HSPM::imposeBlowingVelocity(size_t i, double _blVel)
blVel : double
The blowing velocity to impose. [m/s]
*/
- this->blVel[i] = _blVel;
+ double relax = 1;
+ this->blVel[i] = (1-relax)*this->blVel[i] + relax*_blVel;
}
void HSPM::setdStar(size_t i, double _dStar)
@@ -27,5 +28,6 @@ void HSPM::setdStar(size_t i, double _dStar)
_dStar : double
The dStar value to set. [m]
*/
- this->dStar[i] = _dStar;
+ double relax = 1;
+ this->dStar[i] = (1-relax)*this->dStar[i] + relax*_dStar;
}
\ No newline at end of file
diff --git a/hspm/src/solver.cpp b/hspm/src/solver.cpp
index 9afdfe5ff374dde57649054566e6a7371223fe05..cbfe7226ec5d33a56a28d2afbca0a5effc5ea0fe 100644
--- a/hspm/src/solver.cpp
+++ b/hspm/src/solver.cpp
@@ -27,20 +27,20 @@ void HSPM::solve()
c += blVel;
- // Solve using Eigen's BiCGSTAB solver
- Eigen::BiCGSTAB<Eigen::MatrixXd> solver1;
- solver1.setTolerance(IT_SOLVER_TOLERANCE);
+ // Solve using Eigen's solver
+ Eigen::PartialPivLU<Eigen::MatrixXd> solver1;
+ //solver1.setTolerance(IT_SOLVER_TOLERANCE);
s1 = solver1.compute(A_n).solve(b);
- Eigen::BiCGSTAB<Eigen::MatrixXd> solver2;
- solver2.setTolerance(IT_SOLVER_TOLERANCE);
+ Eigen::PartialPivLU<Eigen::MatrixXd> solver2;
+ //solver2.setTolerance(IT_SOLVER_TOLERANCE);
s2 = solver2.compute(A_n).solve(c);
// Check if the solver converged
- if (solver2.info() != Eigen::Success || solver1.info() != Eigen::Success) {
- std::cout << "Iterative solver failed !" << std::endl;
- exit(-1);
- }
+ //if (solver2.info() != Eigen::Success || solver1.info() != Eigen::Success) {
+ // std::cout << "Iterative solver failed !" << std::endl;
+ // exit(-1);
+ //}
// q = s1*tau + s2
// We need Kutta for tau
@@ -74,6 +74,13 @@ double HSPM::solveOffBodyKutta() {
Eigen::MatrixXd UV_0 = this->getSpeedAtPoint(xc_0, yc_0);
Eigen::MatrixXd UV_N1 = this->getSpeedAtPoint(xc_N1, yc_N1);
+ // Remove the blowing velocity effect
+ // Breaks the 15 deg AoA
+ /*UV_0(0,1) += blVel(0) * sin(theta(0));
+ UV_0(1,1) -= blVel(0) * cos(theta(0));
+ UV_N1(0,1) += blVel(N-1) * sin(theta(N-1));
+ UV_N1(1,1) -= blVel(N-1) * cos(theta(N-1));*/
+
double _a = UV_0(0,0);
double _b = UV_0(0,1);
double _c = UV_0(1,0);
@@ -109,6 +116,18 @@ void HSPM::computeInviscidVelocity()
U(i) = UV(0,0) * tau + UV(0,1);
V(i) = UV(1,0) * tau + UV(1,1);
+
+ // Remove the blowing velocity effect
+ //U(i) += blVel(i) * sin(theta(i));
+ //V(i) -= blVel(i) * cos(theta(i));
+ }
+
+ // In the wake
+ for (size_t i=0; i<N_wake; i++) {
+ Eigen::MatrixXd UV = this->getSpeedAtPoint(x_m_wake(i), y_m_wake(i));
+
+ U_wake(i) = UV(0,0) * tau + UV(0,1);
+ V_wake(i) = UV(1,0) * tau + UV(1,1);
}
}
@@ -133,6 +152,23 @@ void HSPM::computePressureDistribution() {
Cp(i, 1) = y_m(i);
Cp(i, 2) = 0;
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)+1e-10) * sin(theta(i));
+ double y_visc = y_m(i) + (dStar(i)+1e-10) * cos(theta(i));
+ Eigen::MatrixXd UV = this->getSpeedAtPoint(x_visc, y_visc);
+
+ double _U = UV(0,0) * tau + UV(0,1);
+ double _V = UV(1,0) * tau + UV(1,1);
+ double V_mag = sqrt(pow(_U,2) + pow(_V,2));
+
+ // Compute the pressure coefficients
+ Cp(i, 0) = x_m(i);
+ Cp(i, 1) = y_m(i);
+ Cp(i, 2) = 0;
+ Cp(i, 3) = 1 - pow(V_mag / V_inf, 2);
+ */
}
cl = 0;