From f8cdf0c8756124b705f3de81be0ae5a2dd43d70b Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Mon, 17 Dec 2018 01:27:55 -0500
Subject: [PATCH] refactoring thermal conductivity calculation

---
 diffusivityModel/Particle/ParticleI.H         | 26 ++++++-------
 .../diffusivityModel/diffusivityModel.C       | 38 +++++++++++++++++++
 .../diffusivityModel/diffusivityModel.H       |  2 +
 3 files changed, 52 insertions(+), 14 deletions(-)

diff --git a/diffusivityModel/Particle/ParticleI.H b/diffusivityModel/Particle/ParticleI.H
index d5930a9..c8c137c 100644
--- a/diffusivityModel/Particle/ParticleI.H
+++ b/diffusivityModel/Particle/ParticleI.H
@@ -83,15 +83,14 @@ inline Foam::scalar Foam::Particle::Zrot()
 inline Foam::scalar Foam::Particle::zeta(const scalar T, const scalar wellDepth)
     const
 {
-    const scalar a3 = sqrt(pow3(pi));
-    const scalar a1 = a3/2.0;
+    const scalar pi32 = sqrt(pow3(pi));
+
     const scalar a2 = sqr(pi) / 4.0 + 2.0;
 
-    const scalar t1 = sqrt(wellDepth/T);
-    const scalar t2 = wellDepth/T;
-    const scalar t3 = t1 * t2;
+    const scalar rSqrtTstar = 1.0/sqrt(T/wellDepth);
+    const scalar rTstar = wellDepth/T;
 
-    return 1.0 + a1 * t1 + a2 * t2 + a3 * t3;
+    return 1.0 + pi32 * rSqrtTstar * (0.5 + rTstar) + a2 * rTstar;
 }
 
 inline Foam::scalar Foam::Particle::R()
@@ -125,7 +124,7 @@ inline Foam::scalar Foam::Particle::CvVib(const scalar p, const scalar T)
 inline Foam::scalar Foam::Particle::Cv(const scalar p, const scalar T)
     const
 {
-    return thermo_.Cv(p,T) * W();
+    return thermo_.cv(p,T);
 }
 
 
@@ -147,15 +146,14 @@ inline Foam::scalar Foam::Particle::fTrans(const scalar rSc)
 inline Foam::scalar Foam::Particle::lambda(const scalar mu, const scalar D, const scalar Zrot, const GasState &state)
     const
 {
-    const scalar rSc = state.rho() * D / mu;
+    const scalar fVib = state.rho() * D / mu;
 
-    const scalar A = 5./2. - rSc;
-    const scalar B = Zrot + (2./pi) * ((5./3.)*CvRot()/R() + rSc);
-    const scalar AB = (2./pi) * (A/B);
+    const scalar A = 5./2. - fVib;
+    const scalar B = Zrot + (2./pi) * ((5./3.)*CvRot()/R() + fVib);
+    const scalar c1 = (2./pi) * (A/B);
 
-    const scalar fTrans = (5./2.) * (1.0 - AB * CvRot() / CvTrans());
-    const scalar fRot = rSc * (1.0 + AB);
-    const scalar fVib = rSc;
+    const scalar fTrans = (5./2.) * (1.0 - c1 * CvRot() / CvTrans());
+    const scalar fRot = fVib * (1.0 + c1);
 
     return (mu/W()) * (fTrans*CvTrans() + fRot*CvRot() + fVib*CvVib(state.p(), state.T()));
 }
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 433a559..d840384 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -714,6 +714,44 @@ void Foam::diffusivityModel::correct()
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::diffusivityModel::kpure(const scalar T, const scalar p, const UList<scalar> &Y)
+{
+    const speciesTable &species_(thermo_.composition().species());
+
+    scalarSymmetricSquareMatrix Dij(species_.size());
+
+    scalarField localY(Y);
+
+    scalarField Di(species_.size());
+    scalarField Dii(species_.size());
+    scalarField muI(species_.size());
+    scalarField *pkI = new scalarField(species_.size());
+    scalarField &kI = *pkI;
+
+    {
+        const scalar pi = p;
+        const scalar Ti = T;
+
+        GasState state (pi, Ti, localY);
+
+        const scalarField &localX = state.X();
+
+        calculateMuD (muI, Dij, state);
+
+        forAll (Dii, i)
+        {
+            Dii[i] = Dij(i,i);
+        }
+
+        // Pure Thermal conductivities
+        calculateK ( kI, muI, Dii, state);
+
+    }
+
+    return tmp<scalarField>(pkI);
+}
+
+
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 
 void Foam::diffusivityModel::operator=(const diffusivityModel& rhs)
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index 139dd7c..0c4ca8d 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -217,6 +217,8 @@ public:
 
         // Check
 
+        tmp<scalarField> kpure(const scalar T, const scalar p, const UList<scalar> &Y);
+
         // Edit
 
         // Write