From 6e229d946bf1e346eebadc6d74031ceabc643eab Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Sat, 10 Nov 2018 05:00:01 -0500
Subject: [PATCH] Sutherland coefficient calculation using Lennard-Jones
 parameters

---
 ct2foam.C   | 110 ++++++++++++++++++++++++++++++++++++++++++++++++----
 ofFormats.h |   4 +-
 2 files changed, 105 insertions(+), 9 deletions(-)

diff --git a/ct2foam.C b/ct2foam.C
index d1c50e9..9db95cb 100644
--- a/ct2foam.C
+++ b/ct2foam.C
@@ -36,6 +36,8 @@ Description
 
 #include <fstream>
 
+#include <cmath>
+
 #include "cxxopts/cxxopts.hpp"
 
 #include "ofFormats.h"
@@ -54,6 +56,71 @@ std::string stringNASACoefs (doublereal *coef)
     return arr;
 }
 
+void calculateSutherland (Cantera::IdealGasMix *gas_, doublereal As[], doublereal C[])
+{
+    Cantera::Transport *tr_ = Cantera::newTransportMgr("Mix", gas_);
+
+    doublereal mu0[gas_->nSpecies()];
+    doublereal muRatio[gas_->nSpecies()];
+
+    const doublereal T0 = gas_->minTemp();
+    gas_->setState_TP(T0, Cantera::OneAtm);
+    tr_->getSpeciesViscosities(mu0);
+
+    size_t nInterval = 100;
+    doublereal T[nInterval];
+    doublereal TRatioPower[nInterval];
+    doublereal dTemp = (gas_->maxTemp() - gas_->minTemp()) / doublereal(nInterval);
+
+    for (size_t i = 0; i < nInterval; i++)
+    {
+        T[i] = gas_->minTemp() + (i+1) * dTemp;
+    }
+
+    for (size_t i = 0; i < nInterval; i++)
+    {
+        TRatioPower[i] = pow(T[i]/T0, 3./2.);
+    }
+
+    for (size_t j = 0; j < gas_->nSpecies(); j++)
+    {
+        C[j] = 0.0;
+    }
+
+    for (size_t i = 0; i < nInterval; i++)
+    {
+        doublereal Ti = T[i];
+
+        gas_->setTemperature(Ti);
+
+        tr_->getSpeciesViscosities(muRatio);
+        for (size_t j = 0; j < gas_->nSpecies(); j++)
+        {
+            muRatio[j] /= mu0[j];
+        }
+
+        for (size_t j = 0; j < gas_->nSpecies(); j++)
+        {
+            C[j] += (muRatio[j] * Ti - TRatioPower[i] * T0) / (TRatioPower[i] - muRatio[j]);
+        }
+    }
+
+    for (size_t j = 0; j < gas_->nSpecies(); j++)
+    {
+        C[j] /= doublereal(nInterval);
+    }
+
+
+    for (size_t j = 0; j < gas_->nSpecies(); j++)
+    {
+        As[j] = mu0[j] * (T0 + C[j]) / pow(T0, 3./2.);
+    }
+
+    delete tr_;
+
+    return ;
+}
+
 int main(int argc, char *argv[])
 {
 
@@ -105,13 +172,31 @@ int main(int argc, char *argv[])
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Cantera::IdealGasMix gas_(scti, sgas);
-    Cantera::Transport *tr_ = Cantera::newTransportMgr("Mix", &gas_);
-    int nCanteraSp_ = gas_.nSpecies();
 
     fthermo.open(sthermo);
 
     frxn.open(srxn);
 
+    // Thermo Part =============================================================
+
+    // Calculate Sutherland Coefficients
+    doublereal As[gas_.nSpecies()] ;
+    doublereal C[gas_.nSpecies()] ;
+
+    calculateSutherland (&gas_, As, C);
+
+    for (size_t n = 0; n < gas_.nSpecies(); n++)
+    {
+        std::cout << As[n] << std::endl;
+    }
+
+    std::cout << std::endl;
+
+    for (size_t n = 0; n < gas_.nSpecies(); n++)
+    {
+        std::cout << C[n] << std::endl;
+    }
+
     // these constants define the location of coefficient "a6" in the
     // cofficient array c. The c array contains Tmid in the first
     // location, followed by the 7 low-temperature coefficients, then
@@ -135,11 +220,20 @@ int main(int argc, char *argv[])
                     maxTemp,
                     c[0],
                     stringNASACoefs(&c[1]),
-                    stringNASACoefs(&c[8])
+                    stringNASACoefs(&c[8]),
+                    As[n],
+                    C[n]
                     )
                 );
     }
 
+    fthermo.close();
+
+    // End of Thermo Part ======================================================
+
+
+    // Reaction Part
+
     // Write species name list
     frxn<<("species\n");
     frxn<<fmt::format("{}\n", gas_.nSpecies());
@@ -152,7 +246,7 @@ int main(int argc, char *argv[])
     frxn<<(";\n\n");
 
     // Write reaction list
-    frxn<<( "reactions\n{\n");
+    frxn<<( "reactions\n{\n"); // begin "reactions" dictionary
 
     for (size_t k = 0; k < gas_.nReactions(); k++)
     {
@@ -169,7 +263,7 @@ int main(int argc, char *argv[])
         double c[5] = {0};
 
         // reaction name
-        frxn<<( fmt::format( "    un-named-reaction-{}\n", k));
+        frxn<<( fmt::format( "    un-named-reaction-{}\n", k)); // begin a reaction dictionary
         frxn<<( "    {\n" );
 
         switch (gas_.reactionType(k))
@@ -306,10 +400,12 @@ int main(int argc, char *argv[])
                 break;
         }
 
-        frxn<<( "    }\n" );
+        frxn<<( "    }\n" ); // end of a reaction dictionary
     }
 
-    frxn<<( "}\n");
+    frxn<<( "}\n"); // end of "reactions" dictionary
+
+    frxn.close();
 
     return 0;
 }
diff --git a/ofFormats.h b/ofFormats.h
index baf5e1f..fbff9a6 100644
--- a/ofFormats.h
+++ b/ofFormats.h
@@ -22,8 +22,8 @@ std::string thermoFormat = std::string() +
 "    }}\n" +
 "    transport\n" +
 "    {{\n" +
-"        As              1.67212e-06;\n" +
-"        Ts              170.672;\n" +
+"        As              {};\n" + // 8 As, Sutherland Coefficient
+"        Ts              {};\n" + // 9 Ts, Sutherland Temperature
 "    }}\n" +
 "}}\n"
 ;