Updated comment

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;