From aa74d2c6b3096b13c9ef2349f7564cef75335cda Mon Sep 17 00:00:00 2001
From: of <you@example.com>
Date: Sun, 9 Oct 2022 22:28:03 +0800
Subject: [PATCH] update calculation of kappa for single species

---
 src/CanteraMixture/CanteraMixture.C     |   1 +
 src/CanteraMixture/CanteraMixture.H     | 117 ++++++++++++++----------
 src/dfChemistryModel/dfChemistryModel.C |   2 -
 3 files changed, 72 insertions(+), 48 deletions(-)

diff --git a/src/CanteraMixture/CanteraMixture.C b/src/CanteraMixture/CanteraMixture.C
index 706965b3..359b2ad4 100644
--- a/src/CanteraMixture/CanteraMixture.C
+++ b/src/CanteraMixture/CanteraMixture.C
@@ -57,6 +57,7 @@ Foam::CanteraMixture::CanteraMixture
     yTemp_(nSpecies()),
     HaTemp_(nSpecies()),
     CpTemp_(nSpecies()),
+    CvTemp_(nSpecies()),
     muTemp_(nSpecies())
 {
     forAll(Y_, i)
diff --git a/src/CanteraMixture/CanteraMixture.H b/src/CanteraMixture/CanteraMixture.H
index c9160b6c..af0ead42 100644
--- a/src/CanteraMixture/CanteraMixture.H
+++ b/src/CanteraMixture/CanteraMixture.H
@@ -223,71 +223,96 @@ public:
 
 
 private:
-        scalarList HaTemp_;
-        scalarList CpTemp_;
-        scalarList muTemp_;
+
+    scalarList HaTemp_; // J/kmol
+    scalarList CpTemp_; // J/(kmol·k)
+    scalarList CvTemp_; // J/(kmol·k)
+    scalarList muTemp_; // kg/(m·s)
+
+
 public:
 
-        void calcCp(const scalar T, const scalar p)
-        {
-            const scalar RR = constant::physicoChemical::R.value()*1e3; // J/(kmol·k)
+    void calcCp(const scalar T, const scalar p)
+    {
+        const scalar RR = constant::physicoChemical::R.value()*1e3; // J/(kmol·k)
 
-            CanteraGas_->setState_TP(T, p);
+        CanteraGas_->setState_TP(T, p);
 
-            scalarList Cp_R(nSpecies());
-            CanteraGas_->getCp_R(Cp_R.begin());
-            CpTemp_ = Cp_R*RR;
+        scalarList Cp_R(nSpecies());
+        CanteraGas_->getCp_R(Cp_R.begin());
+        CpTemp_ = Cp_R*RR;
+        for(int i=0; i<nSpecies(); ++i)
+        {
+            CvTemp_[i] = CpTemp_[i] - RR;
         }
+    }
 
-        void calcMu(const scalar T, const scalar p)
-        {
-            CanteraGas_->setState_TP(T, p);
+    void calcMu(const scalar T, const scalar p)
+    {
+        CanteraGas_->setState_TP(T, p);
 
-            CanteraTransport_->getSpeciesViscosities(muTemp_.begin());
-        }
+        CanteraTransport_->getSpeciesViscosities(muTemp_.begin());
+    }
 
-        void calcH(const scalar T, const scalar p)
-        {
-            const scalar RT = constant::physicoChemical::R.value()*1e3*T; // J/kmol/K
+    void calcH(const scalar T, const scalar p)
+    {
+        const scalar RT = constant::physicoChemical::R.value()*1e3*T; // J/kmol/K
 
-            CanteraGas_->setState_TP(T, p);
+        CanteraGas_->setState_TP(T, p);
 
-            scalarList Ha_RT(nSpecies());
-            CanteraGas_->getEnthalpy_RT(Ha_RT.begin());
-            HaTemp_ = Ha_RT*RT;
-        }
+        scalarList Ha_RT(nSpecies());
+        CanteraGas_->getEnthalpy_RT(Ha_RT.begin());
+        HaTemp_ = Ha_RT*RT;
+    }
 
-        scalar Cp(label i, scalar p, scalar T) const
-        {
-            return CpTemp_[i]/CanteraGas_->molecularWeight(i);
-        }
+    // J/(kg·K)
+    scalar Cp(label i, scalar p, scalar T) const
+    {
+        return CpTemp_[i]/CanteraGas_->molecularWeight(i);
+    }
 
-        scalar mu(label i, scalar p, scalar T) const
-        {
-            return muTemp_[i];
-        }
+    // J/(kg·K)
+    scalar Cv(label i, scalar p, scalar T) const
+    {
 
-        scalar Ha(label i, scalar p, scalar T) const
-        {
-            return HaTemp_[i]/CanteraGas_->molecularWeight(i);
-        }
+        return CvTemp_[i]/CanteraGas_->molecularWeight(i);
+    }
 
-        scalar Hc(label i) const {return CanteraGas_->Hf298SS(i)/CanteraGas_->molecularWeight(i);} // J/kg
+    // kg/(m·s)
+    scalar mu(label i, scalar p, scalar T) const
+    {
+        return muTemp_[i];
+    }
 
-        scalar Hs(label i, scalar p, scalar T) const {return Ha(i, p, T) - Hc(i);} // J/kg
+    // J/kg
+    scalar Ha(label i, scalar p, scalar T) const
+    {
+        return HaTemp_[i]/CanteraGas_->molecularWeight(i);
+    }
 
-        scalar kappa(label i, scalar p, scalar T) const
-        {
-            // should be kappa of single species
-            // but now lack of access function in Cantera
-            return CanteraTransport_->thermalConductivity();
-        } // W/m/K
+    scalar Hc(label i) const {return CanteraGas_->Hf298SS(i)/CanteraGas_->molecularWeight(i);} // J/kg
+
+    scalar Hs(label i, scalar p, scalar T) const {return Ha(i, p, T) - Hc(i);} // J/kg
+
+    // W/m/K
+    scalar kappa(label i, scalar p, scalar T) const
+    {
+        // should be kappa of single species
+
+        // but now lack of access function in Cantera
+        //return CanteraTransport_->thermalConductivity();
+
+        // use OpenFOAM method to calc kappa
+        scalar Cv_ = Cv(i, p, T); // J/(kg·K)
+        const scalar RR = constant::physicoChemical::R.value()*1e3; // J/(kmol·k)
+        return mu(i, p, T)*Cv_*(1.32 + 1.77*RR/CanteraGas_->molecularWeight(i)/Cv_);
+    } // W/m/K
 
-        scalar Wi(label i) const {return CanteraGas_->molecularWeight(i);} // kg/kmol
+    scalar Wi(label i) const {return CanteraGas_->molecularWeight(i);} // kg/kmol
 
 
-        //- Read dictionary
-        void read(const dictionary&);
+    //- Read dictionary
+    void read(const dictionary&);
 };
 
 
diff --git a/src/dfChemistryModel/dfChemistryModel.C b/src/dfChemistryModel/dfChemistryModel.C
index e80e1a7d..18a10b1e 100644
--- a/src/dfChemistryModel/dfChemistryModel.C
+++ b/src/dfChemistryModel/dfChemistryModel.C
@@ -695,8 +695,6 @@ Foam::dfChemistryModel<ThermoType>::getProblems
 {
     const scalarField& T = T_;
     const scalarField& p = p_;
-    const scalarField& rho = rho_;
-
 
     DynamicList<ChemistryProblem> solved_problems(p.size(), ChemistryProblem(mixture_.nSpecies()));