From 5e1fafa9c748bc50e1f3010d04aca0de7f42ccd5 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Thu, 12 Mar 2020 10:08:42 +0900
Subject: [PATCH] turbulent diffusivity model

---
 .../cantera/transport/TurbulentTransport.h    | 36 +++++++++++++++----
 src/transport/TurbulentTransport.cpp          | 19 +++++-----
 2 files changed, 37 insertions(+), 18 deletions(-)

diff --git a/include/cantera/transport/TurbulentTransport.h b/include/cantera/transport/TurbulentTransport.h
index 89fe7cbf5..69437875c 100644
--- a/include/cantera/transport/TurbulentTransport.h
+++ b/include/cantera/transport/TurbulentTransport.h
@@ -57,10 +57,20 @@ public:
     virtual void getMixDiffCoeffs(double* const d) {
         GasTransport::getMixDiffCoeffs(d);
 
-        double Dm = 0.0;
-        for (size_t k = 0; k < m_nsp; k++) {
-            Dm += d[k] * m_weight[k];
+        doublereal t = m_thermo->temperature();
+        doublereal slope = (-(m_turbmodifier-1.0)/m_profwidth);
+        doublereal trans = slope*(t-m_rdlayerstart) + 1.0;
+
+        double Dm = thermalConductivity() / (m_thermo->density() * m_thermo->cp_mass());
+
+        if (t < m_rdlayerstart - m_profwidth) {
+            Dm *= m_turbmodifier;
+        } else if (t > m_rdlayerstart) {
+            Dm *= 1.0;
+        } else {
+            Dm *= trans;
         }
+
         for (size_t k = 0; k < m_nsp; k++) {
             d[k] = Dm;
         }
@@ -87,10 +97,20 @@ public:
     virtual void getMixDiffCoeffsMass(double* const d){
         GasTransport::getMixDiffCoeffsMass(d);
 
-        double Dm = 0.0;
-        for (size_t k = 0; k < m_nsp; k++) {
-            Dm += d[k] * m_weight[k];
+        doublereal t = m_thermo->temperature();
+        doublereal slope = (-(m_turbmodifier-1.0)/m_profwidth);
+        doublereal trans = slope*(t-m_rdlayerstart) + 1.0;
+
+        double Dm = thermalConductivity() / (m_thermo->density() * m_thermo->cp_mass());
+
+        if (t < m_rdlayerstart - m_profwidth) {
+            Dm *= m_turbmodifier;
+        } else if (t > m_rdlayerstart) {
+            Dm *= 1.0;
+        } else {
+            Dm *= trans;
         }
+
         for (size_t k = 0; k < m_nsp; k++) {
             d[k] = Dm;
         }
@@ -99,7 +119,9 @@ public:
 protected:
 
     //! Internal storage for weights to calculate mean diffusivity
-    vector_fp m_weight;
+    doublereal m_turbmodifier;
+    doublereal m_profwidth;
+    doublereal m_rdlayerstart;
 
 };
 }
diff --git a/src/transport/TurbulentTransport.cpp b/src/transport/TurbulentTransport.cpp
index ace254d50..7b68925bc 100644
--- a/src/transport/TurbulentTransport.cpp
+++ b/src/transport/TurbulentTransport.cpp
@@ -15,7 +15,9 @@ namespace Cantera
 {
 TurbulentTransport::TurbulentTransport() :
     MixTransport(),
-    m_weight(m_nsp)
+    m_turbmodifier(1.0),
+    m_profwidth(1e-5),
+    m_rdlayerstart(300.)
 {
 }
 
@@ -23,18 +25,13 @@ void TurbulentTransport::init(ThermoPhase* thermo, int mode, int log_level)
 {
     MixTransport::init(thermo, mode, log_level);
 
-    m_weight.resize(m_nsp);
-
-    m_weight.assign(m_nsp, 0.0);
-    m_weight[0] = 1.0;
-
     ifstream myfile;
-    myfile.open ("wmean-weights.txt");
-    for (size_t k = 0; k < m_nsp; k++) {
-        myfile >> m_weight[k];
+    myfile.open ("turbulent-coefs.txt");
+
+    myfile >> m_turbmodifier;
+    myfile >> m_profwidth;
+    myfile >> m_rdlayerstart;
 
-        cout << m_thermo->speciesName(k) << m_weight[k] << endl;
-    }
     myfile.close();
 
     // m_reftemp = thermo->temperature();