diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index e3169f49babcdb2935bd6e2bcff0a1c79fa9f1f4..bc83e72f5196fd2dc441306b9ed1f17c52f17b68 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -6,7 +6,7 @@
using namespace tbox;
using namespace dart;
-
+// This is an empty constructor
Closures::Closures()
{
@@ -20,22 +20,22 @@ Closures::~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 delta = (3.15+ H + (1.72/(Hk-1)))*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
+ Hstar = 0.0111*Hks*Hks/dens - 0.0278*Hks*Hks*Hks/dens + 1.528 -0.0002*(Hks*Hk)*(Hks*Hk);
+ Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
+ }
+ 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;
@@ -58,6 +58,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
double Cteq = 0;
double Us = 0;
+
std::vector<double> ClosureVars;
ClosureVars.push_back(Hstar);
@@ -73,35 +74,35 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
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);
+ Hk = std::max(Hk, 1.0005);
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);
+ if(Ret <= 400)
+ {
+ H0 = 4;
+ }
+ else
+ {
+ H0 = 3 + (400/Ret);
+ }
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
+ 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
+ {
+ 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 */
+
+ 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);
@@ -113,35 +114,36 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
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
+ 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
+ else
+ {
+ if (Us > 0.95){
+ Us = 0.98;
+ }
+ n = 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;