Added wake

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++) {

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

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:
     
 };

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));

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

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;