diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index 5d2bae91d01c90df1fd03e8233753424a8794503..b7acfff4d3aa5a69bc1fda57fe636e31ca81a21f 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -30,3 +30,4 @@
#include "wBLRegion.h"
#include "wTimeSolver.h"
#include "wViscFluxes.h"
+#include "wClosures.h"
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index b1eba25b6a021bf5d6c1d3a84b4ea45112ff5420..4e12b24f3c3e837d9e707aa53e39885a56f8d2dc 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -62,6 +62,7 @@ threads="1"
%include "wBLRegion.h"
%include "wTimeSolver.h"
%include "wViscFluxes.h"
+%include "wClosures.h"
%include "wFluid.h"
%include "wAssign.h"
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 72c71064973cf3954e5cdb8e85577529417d6c1b..3a2af7174074169a14ca4270371d7a08eaabbf2d 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -68,6 +68,7 @@ class ViscSolver;
class BLRegion;
class TimeSolver;
class ViscSolver;
+class Closures;
// General
class F0Ps;
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 56115e9bc067e56ae18a2a153a11ce26f0c33449..02cac8eabdff8d7406cc771d674f0a98a8f39717 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -1,4 +1,5 @@
#include "wBLRegion.h"
+#include "wClosures.h"
#include <iostream>
#include <vector>
#include <cmath>
@@ -11,13 +12,16 @@ BLRegion::BLRegion()
{
std::cout << "Coucou de BLRegion\n" << std::endl;
+ closSolver = new Closures();
+
} // end BLRegion
BLRegion::~BLRegion()
-{
- std::cout << "~BLRegion()\n";
+{
+ delete closSolver;
+ std::cout << "~BLRegion()\n";
} // end ~BLRegion
void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
@@ -48,6 +52,12 @@ void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
Regime.resize(nMarkers);
} // End SetMesh
+void BLRegion::SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
+{
+ XxExt = _xxExt;
+ DeltaStarExt = _deltaStarExt;
+}
+
void BLRegion::SetInvBC(std::vector<double> _Ue,
std::vector<double> _Me,
std::vector<double> _Rhoe)
@@ -79,11 +89,58 @@ void BLRegion::SetStagBC(double Re)
U[3] = Ret[0]/(U[0]*Re);
} // End if
- LaminarClosures(0);
+ DeltaStar[0] = U[0] * U[1];
+
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
+
+ Hstar[0] = ClosParam[0];
+ Hstar2[0] = ClosParam[1];
+ Hk[0] = ClosParam[2];
+ Cf[0] = ClosParam[3];
+ Cd[0] = ClosParam[4];
+ Delta[0] = ClosParam[5];
+ Cteq[0] = ClosParam[6];
+ Us[0] = ClosParam[7];
} // End SetStagBC
+void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
+{
+
+ U[0] = UpperCond[0]+LowerCond[0]; /* (Theta_up+Theta_low). */
+ U[1] = (UpperCond[0]*UpperCond[1]+LowerCond[0]*LowerCond[1])/U[0]; /* ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low). */
+ U[2] = 9;
+ U[3] = VtExt[0]; /* Inviscid velocity. */
+ U[4] = (UpperCond[0]*UpperCond[4]+LowerCond[0]*LowerCond[4])/U[0]; /* ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). */
+
+ Ret[0] = U[3]*U[0]*Re; /* Reynolds number based on the momentum thickness. */
+
+ if (Ret[0] < 1){
+ Ret[0] = 1;
+ U[3] = Ret[0] / (U[0]*Re);
+ }
+ DeltaStar[0] = U[0]*U[1]; /* deltaStar */
+
+ /* Laminar closures. */
+
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
+
+ Hstar[0] = ClosParam[0];
+ Hstar2[0] = ClosParam[1];
+ Hk[0] = ClosParam[2];
+ Cf[0] = ClosParam[3];
+ Cd[0] = ClosParam[4];
+ Delta[0] = ClosParam[5];
+ Cteq[0] = ClosParam[6];
+ Us[0] = ClosParam[7];
+
+ for (size_t k=0; k<nMarkers; ++k){
+ Regime[k] = 1;
+ }
+}
+
unsigned int BLRegion::GetnVar(){return nVar;}
size_t BLRegion::GetnMarkers(){return nMarkers;}
+double BLRegion::GetnCrit(){return nCrit;}
double BLRegion::GetLocMarkers(size_t iPoint){return LocMarkers[iPoint];}
double BLRegion::GetGlobMarkers(size_t iPoint){return GlobMarkers[iPoint];}
double BLRegion::GetVtExt(size_t iPoint){return VtExt[iPoint];}
@@ -94,75 +151,11 @@ double BLRegion::GetDeltaStarExt(size_t iPoint){return DeltaStarExt[iPoint];}
void BLRegion::PrintSolution(size_t iPoint)
{
- std::cout << "Point " << iPoint << std::endl;
- std::cout << "Theta " << U[iPoint * nVar + 0] << " \n" << std::endl;
- std::cout << "H " << U[iPoint * nVar + 1] << " \n" << std::endl;
- std::cout << "N " << U[iPoint * nVar + 2] << " \n" << std::endl;
- std::cout << "ue " << U[iPoint * nVar + 3] << " \n" << std::endl;
- std::cout << "Ct " << U[iPoint * nVar + 4] << " \n" << std::endl;
- }
-
-void BLRegion::LaminarClosures(size_t iPoint)
-{
- double theta = U[iPoint*nVar+0];
- double H = U[iPoint*nVar+1];
- double ret = Ret[iPoint];
- double me = Me[iPoint];
- double hk = (H - 0.29*me*me)/(1+0.113*me*me); // Kinematic shape factor
- if(name == "airfoil"){
- hk = std::max(hk, 1.02);
- }
- else{
- hk = std::max(hk, 1.00005);
- }
- hk = std::min(hk,6.0);
- double hstar2 = (0.064/(hk-0.8)+0.251)*me*me; // Density shape parameter
- double delta = std::min((3.15+ H + (1.72/(hk-1)))*theta, 12*theta);
- double hks = hk - 4.35;
- double hstar;
- if(hk <= 4.35){
- double dens = hk + 1;
- hstar = 0.0111*hks*hks/dens - 0.0278*hks*hks*hks/dens + 1.528 -0.0002*(hks*hk)*(hks*hk);
- }
- else if(hk > 4.35){
- hstar = 1.528 + 0.015*hks*hks/hk;
+ std::cout << "Point " << iPoint << std::scientific;
+ std::cout << "Theta " << U[iPoint * nVar + 0] << " \n" << std::scientific;
+ std::cout << "H " << U[iPoint * nVar + 1] << " \n" << std::scientific;
+ std::cout << "N " << U[iPoint * nVar + 2] << " \n" << std::scientific;
+ std::cout << "ue " << U[iPoint * nVar + 3] << " \n" << std::scientific;
+ std::cout << "Ct " << U[iPoint * nVar + 4] << " \n" << std::scientific;
+ std::cout << "Regime " << Regime[iPoint] << " \n" << std::scientific;
}
- hstar = (hstar + 0.028*me*me)/(1+0.014*me*me); // Whitfield's minor additional compressibility correction
-
- /* Cf */
-
- double cf;
- if(hk < 5.5){
- double tmp = (5.5-hk)*(5.5-hk)*(5.5-hk) / (hk+1);
- cf = (-0.07 + 0.0727*tmp)/ret;
- }
- else if(hk >= 5.5){
- double tmp = 1 - 1/(hk-4.5);
- cf = (-0.07 + 0.015*tmp*tmp) /ret;
- }
-
- /* Cd */
-
- double cd;
- if(hk < 4){
- cd = ( 0.00205*pow(4.0-hk, 5.5) + 0.207 ) * (hstar/(2*ret));
- }
- else if(hk >= 4){
- double hkCd = (hk-4)*(hk-4);
- double denCd = 1+0.02*hkCd;
- cd = ( -0.0016*hkCd/denCd + 0.207) * (hstar/(2*ret));
- }
- if (name == "wake"){
- cd = 2*(1.10*(1.0-1.0/hk)*(1.0-1.0/hk)/hk)* (hstar/(2*ret));
- cf = 0;
- }
- Cd[iPoint] = cd;
- Cf[iPoint] = cf;
- Delta[iPoint] = delta;
- Hstar[iPoint] = hstar;
- Hstar2[iPoint] = hstar2;
- Hk[iPoint] = hk;
- Cteq[iPoint] = 0;
- Us[iPoint] = 0;
-}
-
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index 7f0e2ac2e4a8ad6418bbb1f958aba5ebd80bd94c..8d818757e311db4e8d37cb1aaa39ac897819a624 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -2,6 +2,7 @@
#define WBLREGION_H
#include "dart.h"
+#include "wClosures.h"
#include <vector>
#include <string>
@@ -32,10 +33,12 @@ private:
/* Transition related variables */
- double nCrit;
+ double nCrit=9.0;
public:
+ Closures *closSolver;
+
std::string name;
/* Boundary layer variables */
@@ -55,7 +58,7 @@ public:
/* Transition related variables */
- std::vector<double> Regime;
+ std::vector<int> Regime;
unsigned int transMarker;
double xtrF;
@@ -67,11 +70,14 @@ public:
void SetInvBC(std::vector<double> _Ue,
std::vector<double> _Me,
std::vector<double> _Rhoe);
+ void SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt);
void SetStagBC(double Re);
+ void SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond);
unsigned int GetnVar();
size_t GetnMarkers();
+ double GetnCrit();
double GetLocMarkers(size_t iPoint);
double GetGlobMarkers(size_t iPoint);
double GetVtExt(size_t iPoint);
@@ -81,8 +87,6 @@ public:
double GetDeltaStarExt(size_t iPoint);
void PrintSolution(size_t iPoint);
-
- void LaminarClosures(size_t iPoint);
};
} // namespace dart
#endif //WBLREGION_H
\ No newline at end of file
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7eb71544ff0ca171a7e0b743beef6ba58f37ae09
--- /dev/null
+++ b/dart/src/wClosures.cpp
@@ -0,0 +1,159 @@
+#include "wClosures.h"
+#include <iostream>
+#include <math.h>
+#include <string>
+
+using namespace tbox;
+using namespace dart;
+
+
+Closures::Closures()
+{
+ std::cout << "Coucou de Closures\n" << std::endl;
+
+
+} // end Closures
+
+
+Closures::~Closures()
+{
+ std::cout << "~Closures()\n";
+} // end ~Closures
+
+std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name)
+{
+ double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me); // Kinematic shape factor
+ if (name == "upper" || name == "lower"){
+ Hk = std::max(Hk, 1.02);}
+ else{
+ Hk = std::max(Hk, 1.00005);}
+ Hk = std::min(Hk,6.0);
+ double Hstar2 = (0.064/(Hk-0.8)+0.251)*Me*Me; // Density shape parameter
+ double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
+ double Hks = Hk - 4.35;
+ double Hstar;
+ if (Hk <= 4.35){
+ double dens = Hk + 1;
+ Hstar = 0.0111*Hks*Hks/dens - 0.0278*Hks*Hks*Hks/dens + 1.528 -0.0002*(Hks*Hk)*(Hks*Hk);}
+ else if(Hk > 4.35){
+ Hstar = 1.528 + 0.015*Hks*Hks/Hk;}
+ Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
+
+ double tmp;
+ double Cf;
+ double Cd;
+ if (Hk < 5.5){
+ tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1);
+ Cf = (-0.07 + 0.0727*tmp)/Ret;}
+ else if (Hk >= 5.5){
+ tmp = 1 - 1/(Hk-4.5);
+ Cf = (-0.07 + 0.015*tmp*tmp) /Ret;}
+ if (Hk < 4){
+ Cd = ( 0.00205*pow(4.0-Hk, 5.5) + 0.207 ) * (Hstar/(2*Ret));}
+ else if (Hk >= 4){
+ double HkCd = (Hk-4)*(Hk-4);
+ double denCd = 1+0.02*HkCd;
+ Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret));}
+ if (name == "wake"){
+ Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
+ Cf = 0;}
+
+ double Cteq = 0;
+ double Us = 0;
+ std::vector<double> ClosureVars;
+
+ ClosureVars.push_back(Hstar);
+ ClosureVars.push_back(Hstar2);
+ ClosureVars.push_back(Hk);
+ ClosureVars.push_back(Cf);
+ ClosureVars.push_back(Cd);
+ ClosureVars.push_back(delta);
+ ClosureVars.push_back(Cteq);
+ ClosureVars.push_back(Us);
+ return ClosureVars;
+}
+
+std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name)
+{
+ H = std::max(H, 1.0005);
+ // eliminate absurd transients
+ Ct = std::max(std::min(Ct, 0.30), 0.0000001);
+ double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
+ Hk = std::max(std::min(Hk, 10.0), 1.00005);
+ double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
+ //gam = 1.4 - 1
+ //Fc = np.sqrt(1+0.5*gam*Me*Me)
+ double Fc = sqrt(1+0.2*Me*Me);
+
+ double H0;
+ if(Ret <= 400){
+ H0 = 4;}
+ else if(Ret > 400){
+ H0 = 3 + (400/Ret);}
+
+ Ret = std::max(Ret, 200.0);
+
+ double Hstar;
+ if (Hk <= H0){
+ Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;}
+ else if(Hk > H0){
+ Hstar = (Hk-H0)*(Hk-H0)*(0.007*log(Ret)/((Hk-H0+4/log(Ret))*(Hk-H0+4/log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;}
+
+ Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
+ //logRt = log(Ret/Fc);
+ //logRt = std::max(logRt,3);
+ //arg = max(-1.33*Hk, -20)
+ double Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)); // Equivalent normalized wall slip velocity
+ double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
+ double Ctcon = 0.5/(6.7*6.7*0.75);
+
+ double Cf;
+ double Cd;
+ double Hkc;
+ double Cdw;
+ double Cdd;
+ double Cdl;
+ double Cteq;
+ int n;
+ if (name == "wake"){
+ Us = std::min(Us, 0.99995);
+ n = 2; // two wake halves
+ Cf = 0; // no friction inside the wake
+ Hkc = Hk - 1;
+ Cdw = 0; // wall contribution
+ Cdd = (0.995-Us)*Ct*Ct; // turbulent outer layer contribution
+ Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret; // laminar stress contribution to outer layer
+ Cteq = sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
+ }
+ else{
+ if (Us > 0.95){
+ Us = 0.98;
+ }
+ n = 1;
+ //Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
+ Cf = 1/Fc*(0.3*exp(-1.33*Hk)*pow((log10(Ret/Fc)),-1.74-0.31*Hk) + 0.00011*(tanh(4-Hk/0.875) - 1));
+ Hkc = std::max(Hk-1-18/Ret, 0.01);
+ // Dissipation coefficient
+ double Hmin = 1 + 2.1/log(Ret);
+ double Fl = (Hk-1)/(Hmin-1);
+ double Dfac = 0.5+0.5*tanh(Fl);
+ Cdw = 0.5*(Cf*Us) * Dfac;
+ Cdd = (0.995-Us)*Ct*Ct;
+ Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret;
+ Cteq = sqrt(Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
+ }
+ Cd = n*(Cdw+ Cdd + Cdl);
+ // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))*((Hk-1)/(6.432*Hk)));
+
+ std::vector<double> ClosureVars;
+
+ ClosureVars.push_back(Hstar);
+ ClosureVars.push_back(Hstar2);
+ ClosureVars.push_back(Hk);
+ ClosureVars.push_back(Cf);
+ ClosureVars.push_back(Cd);
+ ClosureVars.push_back(delta);
+ ClosureVars.push_back(Cteq);
+ ClosureVars.push_back(Us);
+ return ClosureVars;
+}
diff --git a/dart/src/wClosures.h b/dart/src/wClosures.h
new file mode 100644
index 0000000000000000000000000000000000000000..a51871730f564bc689473f3729f776c0730d1250
--- /dev/null
+++ b/dart/src/wClosures.h
@@ -0,0 +1,20 @@
+#ifndef WCLOSURES_H
+#define WCLOSURES_H
+
+#include "dart.h"
+#include <vector>
+#include <string>
+
+namespace dart
+{
+
+class DART_API Closures
+{
+public:
+ Closures();
+ ~Closures();
+ std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name);
+ std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name);
+};
+} // namespace dart
+#endif //WClosures_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 47c09d27f7fcddd6a5f894669116d849db06217b..dcdd8877b7099e8f9b7e255831572da1fdfb92e8 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -11,17 +11,16 @@ using namespace tbox;
using namespace dart;
-TimeSolver::TimeSolver(double _CFL0, double _Minf, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac)
+TimeSolver::TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac)
{
- std::cout << "Coucou de TimeSolver\n" << std::endl;
CFL0 = _CFL0;
maxIt = _maxIt;
tol = _tol;
itFrozenJac0 = _itFrozenJac;
- Minf = _Minf;
+ Minf = std::max(_Minf, 0.1);
- SysMatrix = new ViscFluxes();
+ SysMatrix = new ViscFluxes(Re);
} // end TimeSolver
@@ -77,6 +76,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
double CFL = CFL0;
unsigned int itFrozenJac = itFrozenJac0;
+
double timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
Matrix<double, 5, 5> IMatrix;
@@ -84,29 +84,53 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
Vector<double, 5> SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
double normSysRes0 = SysRes.norm();
- double normSysRes;
+ double normSysRes = normSysRes0;
Matrix<double, 5, 5> JacMatrix;
Vector<double, 5> slnIncr;
unsigned int nErrors = 0; // Number of errors encountered
unsigned int innerIters = 0; // Inner (non-linear) iterations
+
+ if(bl->name == "upper" && iPoint == 44){
+ std::cout << "IC 44 " << std::endl;
+ bl->PrintSolution(42);
+ }
+
while (innerIters < maxIt){
- if(innerIters % itFrozenJac){
+ if(innerIters % itFrozenJac == 0){
JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
+ if(bl->name == "upper" && iPoint == 44){
+ std::cout << "computeJac at iter " << innerIters << JacMatrix << std::endl;
+ std::cout << "SysRes " << SysRes << std::endl;
+ }
}
-
+
// Ici on peut faire autrement : .asDiagonal()
- slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(SysRes);
+ slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-SysRes);
+
+ if(bl->name == "upper" && iPoint == 44){
+ std::cout << "slnIncr " << slnIncr << std::scientific;}
- for(size_t k=0; k<slnIncr.size(); ++k){
+ for(size_t k=0; k<nVar; ++k){
bl->U[iPoint*nVar+k] += slnIncr(k);
}
SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
normSysRes = (SysRes + slnIncr/timeStep).norm();
+ if(bl->name == "upper" && iPoint == 44){
+ std::cout << "normsysRes/nromsysRes0" << normSysRes/normSysRes0 << " CFL" << CFL << "TimeStep" << timeStep << std::endl;
+ bl->PrintSolution(iPoint);
+ }
+
+ if (normSysRes <= tol * normSysRes0){
+ /* Successfull exit */
+ std::cout << "Point" << iPoint << " Convergence successfull. Niter: " << innerIters << " NormsysRes" << normSysRes << std::endl;
+ return 0;
+ }
+
// CFL adaptation.
CFL = CFL0* pow(normSysRes0/normSysRes, 0.7);
@@ -116,9 +140,9 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
innerIters++;
- } // End while
+ } // End while (innerIters < maxIt)
+ std::cout << "Point" << iPoint << " Max iterations number reached. normSysRes/normSysRes0 " << normSysRes/normSysRes0 << std::endl;
return innerIters;
- return 0;
} // End Integration
double TimeSolver::SetTimeStep(double CFL, double dx, double invVel)
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index b85ac58145e51425189fc0ca6f421c09f2f50cbf..e89628937f3aad81f6d0d9409c130c400653f4b7 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -25,7 +25,7 @@ private:
ViscFluxes *SysMatrix;
public:
- TimeSolver(double _CFL0, double _Minf, unsigned long _maxIt = 1000, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
+ TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 10, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
~TimeSolver();
void InitialCondition(size_t iPoint, BLRegion *bl);
int Integration(size_t iPoint, BLRegion *bl);
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 0374f60d19eea90895555a9e856291bb2f69c4c6..3a66e559520b11cc8d5000a89bd36467c08b1814 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -12,17 +12,16 @@ using namespace dart;
using namespace Eigen;
-ViscFluxes::ViscFluxes()
+ViscFluxes::ViscFluxes(double _Re)
{
- std::cout << "Coucou de ViscFluxes\n" << std::endl;
-
+ Re = _Re;
} // end ViscFluxes
ViscFluxes::~ViscFluxes()
-{
- std::cout << "~ViscFluxes()\n";
+{
+ std::cout << "~ViscFluxes()\n";
} // end ~ViscFluxes
Vector<double, 5> ViscFluxes::ComputeFluxes(size_t iPoint, BLRegion *bl)
@@ -37,6 +36,7 @@ Vector<double, 5> ViscFluxes::ComputeFluxes(size_t iPoint, BLRegion *bl)
Matrix<double, 5, 5> ViscFluxes::ComputeJacobian(size_t iPoint, BLRegion *bl, Vector<double, 5> &sysRes, double eps)
{
+
unsigned int nVar = bl->GetnVar();
Matrix<double, 5, 5> JacMatrix;
std::vector<double> uUp;
@@ -44,8 +44,9 @@ Matrix<double, 5, 5> ViscFluxes::ComputeJacobian(size_t iPoint, BLRegion *bl, Ve
uUp.push_back(bl->U[iPoint*nVar + k]);
} // End for
+ double varSave;
for (size_t iVar=0; iVar<nVar; ++iVar){
- double varSave = uUp[iVar];
+ varSave = uUp[iVar];
uUp[iVar] += eps;
JacMatrix.col(iVar) = (BLlaws(iPoint, bl, uUp) - sysRes) / eps;
@@ -64,22 +65,38 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
double dissipRatio;
if (bl->name == "wake"){
- double dissipRatio = 0.9;
+ dissipRatio = 0.9;
} // End if
else{
- double dissiRatio = 1;
+ dissipRatio = 1;
} // End else
-
- bl->Ret[iPoint] = std::max(u[3] * u[0], 1.0);
+
+ bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0);
bl->DeltaStar[iPoint] = u[0] * u[1];
/* Boundary layer closure relations */
if (bl->Regime[iPoint] == 0){
- LaminarClosures(iPoint, bl);
+ std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[1], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
+ bl->Hstar[iPoint] = lamParam[0];
+ bl->Hstar2[iPoint] = lamParam[1];
+ bl->Hk[iPoint] = lamParam[2];
+ bl->Cf[iPoint] = lamParam[3];
+ bl->Cd[iPoint] = lamParam[4];
+ bl->Delta[iPoint] = lamParam[5];
+ bl->Cteq[iPoint] = lamParam[6];
+ bl->Us[iPoint] = lamParam[7];
} // End if
else if (bl->Regime[iPoint] == 1){
- TurbulentClosures(iPoint, bl);
+ std::vector<double> turbParam = bl->closSolver->TurbulentClosures(u[0], u[1], u[4], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
+ bl->Hstar[iPoint] = turbParam[0];
+ bl->Hstar2[iPoint] = turbParam[1];
+ bl->Hk[iPoint] = turbParam[2];
+ bl->Cf[iPoint] = turbParam[3];
+ bl->Cd[iPoint] = turbParam[4];
+ bl->Delta[iPoint] = turbParam[5];
+ bl->Cteq[iPoint] = turbParam[6];
+ bl->Us[iPoint] = turbParam[7];
} // End else if
else{
@@ -88,30 +105,31 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
/* Gradients computation */
double dx = (bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1));
+ double dxExt = (bl->GetXxExt(iPoint) - bl->GetXxExt(iPoint-1));
double dT_dx = (u[0] - bl->U[(iPoint-1)*nVar + 0])/dx;
double dH_dx = (u[1] - bl->U[(iPoint-1)*nVar + 1])/dx;
double due_dx = (u[3] - bl->U[(iPoint-1)*nVar + 3])/dx;
double dCt_dx = (u[4] - bl->U[(iPoint-1)*nVar + 4])/dx;
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint-1])/dx;
- double dueExt_dx = (bl->GetVtExt(iPoint) - bl->GetVtExt(iPoint))/dx;
- double ddeltaStarExt_dx = (bl->GetDeltaStarExt(iPoint) - bl->GetDeltaStarExt(iPoint))/dx;
+ double dueExt_dx = (bl->GetVtExt(iPoint) - bl->GetVtExt(iPoint-1))/dxExt;
+ double ddeltaStarExt_dx = (bl->GetDeltaStarExt(iPoint) - bl->GetDeltaStarExt(iPoint-1))/dxExt;
double c=4/(M_PI*dx);
- double cExt=4/(M_PI*(bl->GetXxExt(iPoint) - bl->GetXxExt(iPoint-1)));
+ double cExt=4/(M_PI*dxExt);
/* Amplification routine */
double dN_dx;
double Ax;
- if (bl->xtrF!=NULL && bl->Regime[iPoint] == 0){
- double Ax = AmplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
- double dN_dx = (bl->U[iPoint * nVar + 2] - bl->U[(iPoint-1)*nVar + 2])/(bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1));
+ if (bl->xtrF!=-1 && bl->Regime[iPoint] == 0){
+ Ax = AmplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
+ dN_dx = (u[2] - bl->U[(iPoint-1)*nVar + 2])/dx;
} // End if
else{
- double Ax = 0;
- double dN_dx = 0;
+ Ax = 0;
+ dN_dx = 0;
} // End else
double Mea = bl->GetMe(iPoint);
@@ -123,63 +141,92 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
double Cfa = bl->Cf[iPoint];
double Cda = bl->Cd[iPoint];
double deltaa = bl->Delta[iPoint];
+ double Cteqa = bl->Cteq[iPoint];
double Usa = bl->Us[iPoint];
+ /*if (bl->name=="upper" && iPoint == 1){
+ std::cout << "Reta" << Reta << std::endl;
+ std::cout << "Hstara" << Hstara << std::endl;
+ std::cout << "Hstar2a" << Hstar2a << std::endl;
+ std::cout << "Hka" << Hka << std::endl;
+ std::cout << "Cfa" << Cfa << std::endl;
+ std::cout << "Cda" << Cda << std::endl;
+ std::cout << "deltaa" << deltaa << std::endl;
+ std::cout << "Usa" << Usa << std::endl;
+ }*/
+
/* Space part */
- spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0]/u[3]) * due_dx - Cfa;
- spaceVector(1) = u[1] * dHstar_dx + (2*Hstar2a + Hstara * (1-u[1])) * u[0]/u[3] * due_dx - 2*Cda + Hstara * Cfa/2;
+ spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0]/u[3]) * due_dx - Cfa/2;
+ spaceVector(1) = u[0] * dHstar_dx + (2*Hstar2a + Hstara * (1-u[1])) * u[0]/u[3] * due_dx - 2*Cda + Hstara * Cfa/2;
spaceVector(2) = dN_dx - Ax;
spaceVector(3) = due_dx - c*(u[1] * dT_dx + u[0]*dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
if (bl->Regime[iPoint] == 1){
- spaceVector(4) = (2*deltaa/u[4]) * dCt_dx - 5.6*(bl->Cteq[iPoint] - u[4]*dissipRatio)
- - 2*deltaa*(4/(3*u[1]*u[0]) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio)) * ((Hka-1)/(6.7*Hka * dissipRatio))) - 1/u[3] * due_dx);
+ spaceVector(4) = (2*deltaa/u[4]) * dCt_dx - 5.6*(Cteqa - u[4]*dissipRatio)
+ - 2*deltaa*(4/(3*u[1]*u[0]) * (Cfa/2 - (((Hka-1)/(6.7*Hka * dissipRatio)) * ((Hka-1)/(6.7*Hka * dissipRatio)))) - 1/u[3] * due_dx);
} // End if
else{
- spaceVector(4) = 0;
+ spaceVector(4) = 0.0;
}
/* Time part */
timeMatrix(0,0) = u[1]/u[3];
timeMatrix(0,1) = u[0]/u[3];
- timeMatrix(0,2) = 0;
+ timeMatrix(0,2) = 0.0;
timeMatrix(0,3) = u[0]*u[1]/(u[3]*u[3]);
- timeMatrix(0,4) = 0;
+ timeMatrix(0,4) = 0.0;
timeMatrix(1,0) = (1+u[1]*(1-Hstara))/u[3];
- timeMatrix(1,1) = (1-Hstara)*u[1]/u[3];
- timeMatrix(1,2) = 0;
+ timeMatrix(1,1) = (1-Hstara)*u[0]/u[3];
+ timeMatrix(1,2) = 0.0;
timeMatrix(1,3) = (2-Hstara*u[1])*u[0]/(u[3]*u[3]);
- timeMatrix(1,4) = 0;
+ timeMatrix(1,4) = 0.0;
- timeMatrix(2,0) = 0;
- timeMatrix(2,1) = 0;
- timeMatrix(2,2) = 1;
- timeMatrix(2,3) = 0;
- timeMatrix(2,4) = 0;
+ timeMatrix(2,0) = 0.0;
+ timeMatrix(2,1) = 0.0;
+ timeMatrix(2,2) = 1.0;
+ timeMatrix(2,3) = 0.0;
+ timeMatrix(2,4) = 0.0;
timeMatrix(3,0) = -c*u[1];
timeMatrix(3,1) = -c*u[0];
- timeMatrix(3,2) = 0;
- timeMatrix(3,3) = 1;
- timeMatrix(3,4) = 0;
+ timeMatrix(3,2) = 0.0;
+ timeMatrix(3,3) = 1.0;
+ timeMatrix(3,4) = 0.0;
- timeMatrix(4,0) = 0;
- timeMatrix(4,1) = 0;
- timeMatrix(4,2) = 0;
+ timeMatrix(4,0) = 0.0;
+ timeMatrix(4,1) = 0.0;
+ timeMatrix(4,2) = 0.0;
if (bl->Regime[iPoint]==1){
timeMatrix(4,3) = 2*deltaa/(Usa*u[3]*u[3]);
timeMatrix(4,4) = 2*deltaa/(u[3]*u[4]*Usa);
}
else{
- timeMatrix(4,3) = 0;
- timeMatrix(4,4) = 0;
+ timeMatrix(4,3) = 0.0;
+ timeMatrix(4,4) = 1.0;
}
- Vector<double, 5> sysRes = timeMatrix.partialPivLu().solve(spaceVector);
+ /*if(bl->name=="upper" && iPoint==1){
+ std::cout << "timeMatrix" << timeMatrix << std::endl;
+ std::cout << "spaceVector" << spaceVector << std::endl;
+
+ std::cout << "dT_dx" << dT_dx << std::endl;
+ std::cout << "dH_dx" << dH_dx << std::endl;
+ std::cout << "due_dx" << due_dx << std::endl;
+ std::cout << "dCt_dx" << dCt_dx << std::endl;
+ std::cout << "dueExt_dx" << dueExt_dx << std::endl;
+ std::cout << "ddStarExt_dx" << ddeltaStarExt_dx << std::endl;
+ std::cout << "dHstar_dx" << dHstar_dx << std::endl;
+ }*/
+ Vector<double, 5> sysRes;
+ sysRes = timeMatrix.partialPivLu().solve(spaceVector);
+ /*if(bl->name=="upper" && iPoint==1){
+ std::cout << "sysRes" << sysRes << std::endl;
+ }*/
+
return sysRes;
}
@@ -237,163 +284,4 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
}
}
return Ax;
-}
-
-void ViscFluxes::LaminarClosures(size_t iPoint, BLRegion *bl)
-{
- unsigned int nVar = bl->GetnVar();
- double theta = bl->U[iPoint*nVar+0];
- double H = bl->U[iPoint*nVar+1];
- double Ret = bl->Ret[iPoint];
- double Me = bl->GetMe(iPoint);
- double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me); // Kinematic shape factor
- if(bl->name == "airfoil"){
- Hk = std::max(Hk, 1.02);
- }
- else{
- Hk = std::max(Hk, 1.00005);
- }
- Hk = std::min(Hk,6.0);
- double Hstar2 = (0.064/(Hk-0.8)+0.251)*Me*Me; // Density shape parameter
- double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
- double Hks = Hk - 4.35;
- double Hstar;
- if(Hk <= 4.35){
- double dens = Hk + 1;
- Hstar = 0.0111*Hks*Hks/dens - 0.0278*Hks*Hks*Hks/dens + 1.528 -0.0002*(Hks*Hk)*(Hks*Hk);
- }
- else if(Hk > 4.35){
- Hstar = 1.528 + 0.015*Hks*Hks/Hk;
- }
- Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
-
- /* Cf */
-
- double Cf;
- if(Hk < 5.5){
- double tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1);
- Cf = (-0.07 + 0.0727*tmp)/Ret;
- }
- else if(Hk >= 5.5){
- double tmp = 1 - 1/(Hk-4.5);
- Cf = (-0.07 + 0.015*tmp*tmp) /Ret;
- }
-
- /* Cd */
-
- double Cd;
- if(Hk < 4){
- Cd = ( 0.00205*pow(4.0-Hk, 5.5) + 0.207 ) * (Hstar/(2*Ret));
- }
- else if(Hk >= 4){
- double HkCd = (Hk-4)*(Hk-4);
- double denCd = 1+0.02*HkCd;
- Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret));
- }
- if (bl->name == "wake"){
- Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
- Cf = 0;
- }
- bl->Cd[iPoint] = Cd;
- bl->Cf[iPoint] = Cf;
- bl->Delta[iPoint] = delta;
- bl->Hstar[iPoint] = Hstar;
- bl->Hstar2[iPoint] = Hstar2;
- bl->Hk[iPoint] = Hk;
- bl->Cteq[iPoint] = 0;
- bl->Us[iPoint] = 0;
-}
-
-void ViscFluxes::TurbulentClosures(size_t iPoint, BLRegion *bl)
-{
- unsigned int nVar = bl->GetnVar();
- double theta = bl->U[iPoint*nVar+0];
- double H = bl->U[iPoint*nVar+1];
- double Ct = bl->U[iPoint*nVar+4];
- double Ret = bl->Ret[iPoint];
- double Me = bl->GetMe(iPoint);
-
- H = std::max(H, 1.0005);
- /* Eliminate absurd transients */
- Ct = std::max(std::min(Ct, 0.30), 0.0000001);
- double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
- Hk = std::max(std::min(Hk, 10.0), 1.00005);
- double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
- //gam = 1.4 - 1
- //Fc = np.sqrt(1+0.5*gam*Me**2)
- double Fc = sqrt(1+0.2*Me*Me);
-
- double H0;
- if(Ret <= 400){
- H0 = 4;
- }
- else if(Ret > 400){
- H0 = 3 + (400/Ret);
- }
-
- Ret = std::max(Ret, 200.0);
- double Hstar;
- if (Hk <= H0){
- Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;
- }
- else if (Hk > H0){
- Hstar = (Hk-H0)*(Hk-H0)*(0.007*log(Ret)/((Hk-H0+4/log(Ret))*(Hk-H0+4/log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
- }
-
- Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
- //logRt = np.log(Ret/Fc)
- //logRt = max(logRt,3)
- //arg = max(-1.33*Hk, -20)
- double Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)); // Equivalent normalized wall slip velocity
- double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
- double Ctcon = 0.5/(6.7*6.7*0.75);
-
- int n;
- double Cf;
- double Hkc;
- double Cdw;
- double Cdd;
- double Cdl;
- double Cteq;
- double Hmin;
- double Fl;
- double Dfac;
-
- if (bl->name == "wake"){
- Us = std::min(Us, 0.99995);
- n = 2; // two wake halves
- Cf = 0; // no friction inside the wake
- Hkc = Hk - 1;
- Cdw = 0; // wall contribution
- Cdd = (0.995-Us)*Ct*Ct; // turbulent outer layer contribution
- Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret; // laminar stress contribution to outer layer
- Cteq = sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
- }
- else{
- if(Us > 0.95){
- Us = 0.98;
- } // End if
- n = 1;
- // Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
- Cf = 1/Fc*(0.3*exp(-1.33*Hk)*pow(log10(Ret/Fc),-1.74-0.31*Hk) + 0.00011*(tanh(4-Hk/0.875) - 1));
- Hkc = std::max(Hk-1-18/Ret, 0.01);
- // Dissipation coefficient
- Hmin = 1 + 2.1/log(Ret);
- Fl = (Hk-1)/(Hmin-1);
- Dfac = 0.5+0.5*tanh(Fl);
- Cdw = 0.5*(Cf*Us) * Dfac;
- Cdd = (0.995-Us)*Ct*Ct;
- Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret;
- Cteq = sqrt(Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
- } // End else
- double Cd = n*(Cdw+ Cdd + Cdl);
- // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk)*((Hk-1)/(6.432*Hk)));
- bl->Cd[iPoint] = Cd;
- bl->Cf[iPoint] = Cf;
- bl->Delta[iPoint] = delta;
- bl->Hstar[iPoint] = Hstar;
- bl->Hstar2[iPoint] = Hstar2;
- bl->Hk[iPoint] = Hk;
- bl->Cteq[iPoint] = Cteq;
- bl->Us[iPoint] = Us;
}
\ No newline at end of file
diff --git a/dart/src/wViscFluxes.h b/dart/src/wViscFluxes.h
index 8e9ac8ef9fe59096aa9b93e8c3cd5d4521f83b9a..f79f87ea3fa05c2d605744f44484a5f3c36c1eaf 100644
--- a/dart/src/wViscFluxes.h
+++ b/dart/src/wViscFluxes.h
@@ -12,16 +12,17 @@ namespace dart
class DART_API ViscFluxes
{
+private:
+ double Re;
public:
- ViscFluxes();
+ ViscFluxes(double _Re);
~ViscFluxes();
Eigen::Vector<double, 5> ComputeFluxes(size_t iPoint, BLRegion *bl);
Eigen::Matrix<double, 5, 5> ComputeJacobian(size_t iPoint, BLRegion *bl, Eigen::Vector<double, 5> &sysRes, double eps);
private:
Eigen::Vector<double, 5> BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u);
double AmplificationRoutine(double hk, double ret, double theta);
- void LaminarClosures(size_t iPoint, BLRegion *bl);
- void TurbulentClosures(size_t iPoint, BLRegion *bl);
+
};
} // namespace dart
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 0ba78fd28f634469afb37ffa462b03d446e81151..09db6df85b8de976eb8bc4f9ac4187bb4858052a 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -19,9 +19,6 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
/* User input numerical parameters */
CFL0 = _CFL0;
- /* Parameters */
- double cd = 0.0;
-
std::cout << "--- Viscous solver setup ---\n"
<< "Reynolds number: " << Re
<< "Freestream Mach number: " << Minf
@@ -32,17 +29,13 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
//std::vector<dart::BLRegion> bl;
- std::cout << "--- Jusqu'ici ok ---\n" << std::endl;
-
for (size_t i=0; i<3; ++i)
{
- std::cout << "yop" << std::endl;
BLRegion *region = new BLRegion();
- std::cout << "yop2" << std::endl;
bl.push_back(region);
} // end for
- tSolver = new TimeSolver(CFL0, Minf);
+ tSolver = new TimeSolver(CFL0, Minf, Re);
bl[0]->name = "upper";
bl[1]->name = "lower";
@@ -71,18 +64,30 @@ void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector
bl[region]->SetMesh(locM, globM);
}
-void ViscSolver::SetInvBC(size_t region, std::vector<double> _Ue,
+void ViscSolver::SendInvBC(size_t region, std::vector<double> _Ue,
std::vector<double> _Me,
std::vector<double> _Rhoe)
{
bl[region]->SetInvBC(_Ue, _Me, _Rhoe);
}
+void ViscSolver::SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
+{
+ bl[region]->SetExtVariables(_xxExt, _deltaStarExt);
+}
-int ViscSolver::Run(){
+
+int ViscSolver::Run(unsigned int couplIter){
std::cout << "Entered run"
<< std::endl;
+ /* Prepare time solver */
+
+ if (couplIter == 0){
+ tSolver->SetCFL0(1);
+ tSolver->SetitFrozenJac(1);
+ }
+
/* Set boundary condition */
bl[0]->SetStagBC(Re);
@@ -90,21 +95,188 @@ int ViscSolver::Run(){
/* Loop over the regions */
+ bool lockTrans;
+
for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion){
+ /* Reset regime */
+ if (iRegion<2){
+ for (size_t k=0; k<bl[iRegion]->GetnMarkers(); k++){
+ bl[iRegion]->Regime[k]=0;
+ }
+ }
+ else{
+ for (size_t k=0; k<bl[iRegion]->GetnMarkers(); k++){
+ bl[iRegion]->Regime[k]=1;
+ }
+ }
+
+ lockTrans = false;
+
+ std::cout << "Starting " << bl[iRegion]->name << std::endl;
+
+ if (iRegion == 2){
+ /* Impose wake boundary condition */
+ std::vector<double> UpperCond;
+ std::vector<double> LowerCond;
+ UpperCond = std::vector<double>(bl[0]->U.end() - bl[0]->GetnVar(), bl[0]->U.end()); /* Base std::vector constructor to slice */
+ LowerCond = std::vector<double>(bl[1]->U.end() - bl[1]->GetnVar(), bl[1]->U.end()); /* Base std::vector constructor to slice */
+ bl[iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
+
+ lockTrans = true;
+ }
+
/* Loop over points */
for (size_t iPoint = 1; iPoint < bl[iRegion]->GetnMarkers(); ++iPoint){
+ /* Initialize solution at point */
+
tSolver->InitialCondition(iPoint, bl[iRegion]);
- std::cout << "Point " << iPoint << " H " << bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 3] << std::endl;
+ /* Time integration */
+
+ tSolver->Integration(iPoint, bl[iRegion]);
+
+ /* Transition */
+
+ if (!lockTrans){
+ CheckWaves(iPoint, bl[iRegion]);
+
+ // Free transition.
+ if (bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 2] >= bl[iRegion]->GetnCrit()){
+ AverageTransition(iPoint, bl[iRegion], 0);
+ std::cout << "Free transition x/c = " << bl[iRegion]->xtr << std::endl;
+ lockTrans = true;
+ } // End if
+ } // End if (!lockTrans).
+
+ } // End for iPoints
+ } // End for iRegions
+
+ ComputeDrag();
+ return 0; // Successfull exit
+}
+
+void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
+{
+ /* Determine if the amplification waves start growing or not. */
- //tSolver->Integration(iPoint, bl[iRegion]);
+ double logRet_crit = 2.492*pow(1/(bl->Hk[iPoint]-1), 0.43)
+ + 0.7*(tanh(14*(1/(bl->Hk[iPoint]-1))-9.24)+1.0);
+ if (bl->Ret[iPoint] > 0){ // Avoid log of negative number.
+ if (log10(bl->Ret[iPoint]) <= logRet_crit - 0.08){
+ bl->U[iPoint*bl->GetnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
}
+ }
+}
-
+void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
+{
+ /* Averages solution on transition marker */
+
+ unsigned int nVar = bl->GetnVar();
+
+ // Save laminar solution.
+ std::vector<double> lamSol;
+ for (size_t k=0; k<nVar; ++k){
+ lamSol.push_back(bl->U[iPoint*nVar + k]);
}
- return 0; // Successfull exit
-};
\ No newline at end of file
+
+ // Set up turbulent portion boundary condition.
+ bl->transMarker = iPoint; /* Save transition marker */
+
+ /* Loop over remaining points in the region */
+ for (size_t k=0; k<bl->GetnMarkers() - iPoint; ++k){
+ bl->Regime[iPoint+k] = 1; // Set Turbulent regime.
+ }
+
+ /* Compute transition location */
+ bl->xtr = (bl->GetnCrit()-bl->U[(iPoint-1)*nVar + 2])*((bl->GetGlobMarkers(iPoint)-bl->GetGlobMarkers(iPoint-1))/(bl->U[iPoint*nVar + 2]-bl->U[(iPoint-1)*nVar+2]))+bl->GetGlobMarkers(iPoint-1);
+
+ /* Percentage of the subsection corresponding to a turbulent flow */
+ double avTurb = (bl->GetGlobMarkers(iPoint)-bl->xtr)/(bl->GetGlobMarkers(iPoint)-bl->GetGlobMarkers(iPoint-1));
+ double avLam = 1-avTurb;
+
+ /* Impose boundary condition */
+ std::vector<double> turbParam1 = bl->closSolver->TurbulentClosures(bl->U[(iPoint-1)*nVar + 0], bl->U[(iPoint-1)*nVar + 1], 0.03, bl->Ret[iPoint-1], bl->GetMe(iPoint-1), bl->name);
+ bl->Cteq[iPoint-1] = turbParam1[6];
+ bl->U[(iPoint-1)*nVar + 4] = 0.7 * bl->Cteq[iPoint-1];
+
+ /* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched line 183) */
+ /* These initial guess values do not influence the converged solution */
+ bl->U[iPoint*nVar + 4] = bl->U[(iPoint-1)*nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
+ bl->U[iPoint*nVar + 1] = 1.515; // Because we expect a significant drop of the shape factor.
+
+ /* Solve point in turbulent configuration */
+ tSolver->Integration(iPoint, bl);
+
+ /* Average both solutions (Now solution @ iPoint is the turbulent solution) */
+ bl->U[iPoint*nVar + 0] = avLam * lamSol[0] + avTurb * bl->U[(iPoint)*nVar + 0]; /* Theta */
+ bl->U[iPoint*nVar + 1] = avLam * lamSol[1] + avTurb * bl->U[(iPoint)*nVar + 1]; /* H */
+ bl->U[iPoint*nVar + 2] = 9; /* N */
+ bl->U[iPoint*nVar + 3] = avLam * lamSol[3] + avTurb * bl->U[(iPoint)*nVar + 3]; /* ue */
+ bl->U[iPoint*nVar + 4] = avTurb * bl->U[(iPoint)*nVar + 4]; /* Ct */
+
+ /* Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). */
+ std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[(iPoint-1)*nVar + 0], bl->U[(iPoint-1)*nVar + 1], bl->Ret[iPoint-1], bl->GetMe(iPoint-1), bl->name);
+ bl->Hstar[iPoint-1] = lamParam[0];
+ bl->Hstar2[iPoint-1] = lamParam[1];
+ bl->Hk[iPoint-1] = lamParam[2];
+ bl->Cf[iPoint-1] = lamParam[3];
+ bl->Cd[iPoint-1] = lamParam[4];
+ bl->Delta[iPoint-1] = lamParam[5];
+ bl->Cteq[iPoint-1] = lamParam[6];
+ bl->Us[iPoint-1] = lamParam[7];
+
+ std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[(iPoint)*nVar + 0], bl->U[(iPoint)*nVar + 1], bl->U[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
+ bl->Hstar[iPoint] = turbParam[0];
+ bl->Hstar2[iPoint] = turbParam[1];
+ bl->Hk[iPoint] = turbParam[2];
+ bl->Cf[iPoint] = turbParam[3];
+ bl->Cd[iPoint] = turbParam[4];
+ bl->Delta[iPoint] = turbParam[5];
+ bl->Cteq[iPoint] = turbParam[6];
+ bl->Us[iPoint] = turbParam[7];
+}
+
+
+void ViscSolver::ComputeDrag()
+{
+
+ unsigned int nVar = bl[0]->GetnVar();
+ /* Normalize friction and dissipation coefficients. */
+
+ /* std::vector<double> frictioncoeff_upper;
+ for (size_t k=0; k<bl[0]->GetnMarkers(); ++k){
+ frictioncoeff_upper.push_back(bl[0]->U[k*nVar + 3] * bl[0]->Cf[k]);
+ }
+ std::vector<double> frictioncoeff_lower;
+ for (size_t k=0; k<bl[1]->GetnMarkers(); ++k){
+ frictioncoeff_upper.push_back(bl[1]->U[k*nVar + 3] * bl[1]->Cf[k]);
+ }
+
+ /* Compute friction drag coefficient. */
+
+ /* double Cdf_upper = 0.0;
+ for (size_t k=1; k<bl[0]->GetnMarkers(); ++k){
+ Cdf_upper += (bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k-1))*frictioncoeff_upper[k-1];
+ }
+ double Cdf_lower = 0.0;
+ for (size_t k=1; k<bl[1]->GetnMarkers(); ++k){
+ Cdf_upper += (bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k-1))*frictioncoeff_lower[k-1];
+ }
+
+ */
+
+ // Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
+
+ /* Compute total drag coefficient. */
+
+ Cdt = 2 * bl[2]->U[bl[2]->U.back()-5] * pow(bl[2]->U[bl[2]->U.back()-2],(bl[2]->U[bl[2]->U.back()-4]+5)/2);
+
+ /* Compute pressure drag coefficient. */
+
+ //Cdp = Cdt - Cdf;
+}
\ No newline at end of file
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index 5a369f480d4b07620b19963dbe2fd9aeb6bc8290..0338f9fe7dd866f4bfbbf5d1952e096de3680c2b 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -10,7 +10,9 @@ namespace dart
class DART_API ViscSolver
{
public:
- double cd;
+ double Cdt;
+ double Cdf;
+ double Cdp;
private:
double Re;
@@ -26,12 +28,18 @@ public:
~ViscSolver();
void InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM);
- void SetInvBC(size_t region, std::vector<double> _Ue,
+ void SendInvBC(size_t region, std::vector<double> _Ue,
std::vector<double> _Me,
std::vector<double> _Rhoe);
+ void SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt);
- int Run();
+ int Run(unsigned int couplIter);
+
+private:
+ void CheckWaves(size_t iPoint, BLRegion *bl);
+ void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
+ void ComputeDrag();
};
} // namespace dart
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 674660638eb822e759637b6cd0f6c4ebe166d8ad..80f4942b5b2b87262a6c07876d3b9580e00b52a6 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -73,7 +73,7 @@ def main():
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 5
+ CFL0 = 1
# ========================================================================================
@@ -86,7 +86,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
#pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index 8391d1f52c94c7d5bca46afb408361614aaf1d98..6543ced2e408a8b12f6b8e6e38bb183e45ee1d2b 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -68,7 +68,7 @@ def main():
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 5
+ CFL0 = 1
# ========================================================================================
@@ -112,7 +112,7 @@ def main():
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' Re M alpha Cl Cd Cdp Cdf Cm')
- print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
# display timers
tms['total'].stop()
@@ -124,10 +124,10 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
- tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+ tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
- val=validation.Validation(case)
- val.main(isolver.Cl,vsolver.wData)
+ #val=validation.Validation(case)
+ #val.main(isolver.Cl,vsolver.wData)
# check results
@@ -135,7 +135,7 @@ def main():
tests = CTests()
if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
- tests.add(CTest('Cd', vsolver.Cd, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
+ tests.add(CTest('Cd', vsolver.Cdt, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 8a575f0732bb7db28b8293c55acd9b8d1406e996..75253b5598709e1cdde173133219624561c8348e 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -32,7 +32,7 @@ class Coupler:
self.iSolver=iSolver
self.vSolver=vSolver
- self.maxCouplIter = 2
+ self.maxCouplIter = 10
self.couplTol = 1e-4
pass
@@ -75,28 +75,37 @@ class Coupler:
self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.InitMeshes(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.InitMeshes(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['globCoord'][fMeshDict['bNodes'][2]:])
+
+ self.vSolver.SendExtVariables(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
+ self.vSolver.SendExtVariables(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ self.vSolver.SendExtVariables(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
print('fro mpy : Mesh ok')
- self.vSolver.SetInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
- self.vSolver.SetInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
- self.vSolver.SetInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
+ else:
+ self.vSolver.SendExtVariables(0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
+ self.vSolver.SendExtVariables(1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ self.vSolver.SendExtVariables(2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
+
+ self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
+ self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
print('from py: InvBC ok')
# Run viscous solver
- exitCode = self.vSolver.Run()
+ exitCode = self.vSolver.Run(couplIter)
print('From python exitcode: ', exitCode)
- error = (self.vSolver.cd - cdPrev)/self.vSolver.cd if self.vSolver.cd != 0 else 1
+ error = (self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt if self.vSolver.Cdt != 0 else 1
if error <= self.couplTol:
print('')
print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
- .format(couplIter, self.vsolver.cd, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
+ .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
print('')
return 0
print('')
print(ccolors.ANSI_BLUE, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} > {:<2.3f}'
- .format(couplIter, self.vSolver.cd, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
+ .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
print('')
# Set inviscid boundary condition
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index cb6d02c925a2736334e5110b99580da1b919a2a2..176a67873cd83849814e28c1c479679880c53b2e 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -1,19 +1,3 @@
-################################################################################
-# #
-# Flow viscous module file #
-# Constructor : Merges and sorts data #
-# of the 3 groups (upper, lower sides and wake) #
-# #
-# Author: Paul Dechamps #
-# Université de Liège #
-# Date: 2021.09.10 #
-# #
-################################################################################
-
-
-# ------------------------------------------------------------------------------
-# Imports
-# ------------------------------------------------------------------------------
from matplotlib import pyplot as plt
import numpy as np
@@ -42,38 +26,38 @@ class GroupConstructor():
if mapMesh:
# Save initial mesh.
- cMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
- cMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0])]
+ cMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ cMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
- cxxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][:,8],dataIn[1][:,8]))
+ cxxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
# Create fine mesh.
fMeshUp = self.createFineMesh(dataIn[0][0][:,0], nFactor)
fMeshLw = self.createFineMesh(dataIn[0][1][:,0], nFactor)
fMeshWk = self.createFineMesh(dataIn[1][:,0] , nFactor)
- fMesh = np.concatenate((fMeshUp, fMeshLw, fMeshWk))
- fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw)]
+ fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
+ fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
fxx_up = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
fxx_lw = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
fxx_wk = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
- fxx = np.concatenate((fxx_up, fxx_lw, fxx_wk))
+ fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
fvtExt_up = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
fvtExt_lw = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
fvtExt_wk = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
- fvtExt = np.concatenate((fvtExt_up, fvtExt_lw, fvtExt_wk))
+ fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
fMe_up = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
fMe_lw = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
fMe_wk = self.interpolateToFineMesh(dataIn[1][:,5] , fMeshWk, nFactor)
- fMe = np.concatenate((fMe_up, fMe_lw, fMe_wk))
+ fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
frho_up = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
frho_lw = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
frho_wk = self.interpolateToFineMesh(dataIn[1][:,6] , fMeshWk, nFactor)
- frho = np.concatenate((frho_up, frho_lw, frho_wk))
+ frho = np.concatenate((frho_up, frho_lw[1:], frho_wk))
fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
@@ -83,17 +67,17 @@ class GroupConstructor():
fxxExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
fxxExt_wk = self.interpolateToFineMesh(dataIn[1][:,8] , fMeshWk, nFactor)
- fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw, fdStarExt_wk))
- fxxext = np.concatenate((fxxExt_up, fxxExt_lw, fxxExt_wk))
+ fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
+ fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
- fdv = np.concatenate((dv_up, dv_lw, dv_wk))
+ fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
- cMe = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][:,5], dataIn[1][:,5]))
- cRho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][:,6], dataIn[1][:,6]))
- cdStarExt = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][:,7], dataIn[1][:,7]))
+ cMe = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+ cRho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+ cdStarExt = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
@@ -113,18 +97,18 @@ class GroupConstructor():
else:
Mesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
- MeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+ MeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0])]
- xx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
- vtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+ xx = np.concatenate((xx_up, xx_lw, xx_wk))
+ vtExt = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk))
alpha = np.concatenate((alpha_up, alpha_lw, alpha_wk))
- dv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+ dv = np.concatenate((dv_up, dv_lw, dv_wk))
- Me = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
- rho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
- dStarext = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
- xxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8], dataIn[1][:,8]))
+ Me = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][:,5], dataIn[1][:,5]))
+ rho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][:,6], dataIn[1][:,6]))
+ dStarext = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][:,7], dataIn[1][:,7]))
+ xxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][:,8], dataIn[1][:,8]))
fMeshDict = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
@@ -187,4 +171,4 @@ class GroupConstructor():
fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
fMesh = np.append(fMesh, cMesh[-1])
- return fMesh
\ No newline at end of file
+ return fMesh
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
index 423b169411de9d3a6fdbeb0fcdd95a82e8737d03..4c97ee0622c8ca7bf34d4f6b135b575d1252651a 100644
--- a/dart/viscous/Solvers/DIntegration.py
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -108,9 +108,9 @@ class Integration:
nErrors = 0 # Count for the number of errors identified.
innerIters = 0 # Number of inner iterations to solver the non-linear system.
while innerIters <= self.__maxIt:
-
try:
-
+ if iMarker == 2:
+ quit()
# Jacobian computation.
if innerIters % itFrozenJac == 0:
@@ -129,6 +129,7 @@ class Integration:
sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
normSysRes = np.linalg.norm(sysRes + slnIncr/timeStep)
+ print('CFL', CFL, 'timeStep', timeStep)
convPlot.append(normSysRes/normSysRes0)
# Check convergence.
@@ -142,6 +143,7 @@ class Integration:
quit()
except Exception:
+ print('it ', innerIters, ' Entered Exception')
nErrors += 1
diff --git a/dart/viscous/Solvers/DSysMatrix.py b/dart/viscous/Solvers/DSysMatrix.py
index 1da7e742c2434f4bf2a4c3a851a04e37f78b26e9..f19b0819a36ae160d00d521a4aa81279bf234847 100644
--- a/dart/viscous/Solvers/DSysMatrix.py
+++ b/dart/viscous/Solvers/DSysMatrix.py
@@ -41,6 +41,8 @@ class flowEquations:
- linSysRes (np.ndarray): Residuals of the linear system at the considered point.
"""
+ #print('IC', x)
+
timeMatrix = np.empty((dataCont.nVar, dataCont.nVar))
spaceMatrix = np.empty(dataCont.nVar)
@@ -71,6 +73,15 @@ class flowEquations:
Mea = dataCont.extPar[iMarker*dataCont.nExtPar+0]
Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
+ """if iMarker == 1:
+ print("Reta",dataCont.blPar[iMarker*dataCont.nBlPar+0])
+ print("Hstara",Hstara)
+ print("Hstar2a",Hstar2a)
+ print("Hka",Hka)
+ print("Cfa",Cfa)
+ print("Cda",Cda)
+ print("deltaa",deltaa)
+ print("Usa",Usa)"""
# Amplification factor for e^n method.
if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
@@ -145,6 +156,21 @@ class flowEquations:
timeMatrix[4] = [0, 0, 0, 0, 1]
sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
+
+ """print('dTheta_dx',dTheta_dx)
+ print('dH_dx',dH_dx)
+ print('due_dx',due_dx)
+ print('dCt_dx',dCt_dx)
+ print('dueExt_dx',dueExt_dx)
+ print('ddstarExt_dx',ddeltaStarExt_dx)
+ print('dHstar_dx', dHstar_dx)
+ print('2ndTerm', (2+Ha-Mea**2) * (Ta/uea) * due_dx)
+ print('third term', Cfa/2)
+ print('cteq', Cteqa)
+
+ print('timeMatrix', timeMatrix)
+ print('spaceVector', spaceMatrix)
+ print('sysRes', sysRes)"""
"""sysRes = np.empty(5)
if dataCont.flowRegime[iMarker] == 0: # Laminar case