From 6167a724e5b65e6ad52e5b9343f57b35c66a7ab1 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Fri, 11 May 2012 15:27:56 +0000
Subject: [PATCH] Added functions for computing alternate mixture-averaged
 diffusion coefficients

---
 include/cantera/transport/GasTransport.h | 23 ++++++++
 src/transport/GasTransport.cpp           | 67 +++++++++++++++++++++++-
 2 files changed, 88 insertions(+), 2 deletions(-)

diff --git a/include/cantera/transport/GasTransport.h b/include/cantera/transport/GasTransport.h
index 6721a3b2b..d83863d5c 100644
--- a/include/cantera/transport/GasTransport.h
+++ b/include/cantera/transport/GasTransport.h
@@ -76,6 +76,29 @@ public:
      */
     virtual void getMixDiffCoeffs(doublereal* const d);
 
+    //! Returns the mixture-averaged diffusion coefficients [m^2/s].
+    //! These are the coefficients for calculating the molar diffusive fluxes
+    //! from the species mole fraction gradients, computed according to
+    //! Eq. 12.176 in "Chemically Reacting Flow":
+    //!
+    //! \f[  D_{km}^* = \frac{1-X_k}{\Sum_{j \ne k}^K X_j/\mathcal{D}_{kj}} \f]
+    //!
+    //! @param[out] d vector of mixture-averaged diffusion coefficients for
+    //!     each species, length m_nsp.
+    virtual void getMixDiffCoeffsMole(doublereal* const d);
+
+    //! Returns the mixture-averaged diffusion coefficients [m^2/s].
+    //! These are the coefficients for calculating the diffusive mass fluxes
+    //! from the species mass fraction gradients, computed according to
+    //! Eq. 12.178 in "Chemically Reacting Flow":
+    //!
+    //! \f[  \frac{1}{D_{km}} = \Sum_{j \ne k}^K \frac{X_j}{\mathcal{D}_{kj}} +
+    //!     \frac{X_k}{1-Y_k} \Sum_{j \ne k}^K \frac{Y_j}{\mathcal{D}_{kj}} \f]
+    //!
+    //! @param[out] d vector of mixture-averaged diffusion coefficients for
+    //!     each species, length m_nsp.
+    virtual void getMixDiffCoeffsMass(doublereal* const d);
+
 protected:
     GasTransport(ThermoPhase* thermo=0);
 
diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp
index c1d6dee14..371c05f1e 100644
--- a/src/transport/GasTransport.cpp
+++ b/src/transport/GasTransport.cpp
@@ -273,7 +273,7 @@ void GasTransport::getMixDiffCoeffs(doublereal* const d)
     }
 
     doublereal mmw = m_thermo->meanMolecularWeight();
-    doublereal sumxw = 0.0, sum2;
+    doublereal sumxw = 0.0;
     doublereal p = m_thermo->pressure();
     if (m_nsp == 1) {
         d[0] = m_bdiff(0,0) / p;
@@ -282,7 +282,7 @@ void GasTransport::getMixDiffCoeffs(doublereal* const d)
             sumxw += m_molefracs[k] * m_mw[k];
         }
         for (size_t k = 0; k < m_nsp; k++) {
-            sum2 = 0.0;
+            double sum2 = 0.0;
             for (size_t j = 0; j < m_nsp; j++) {
                 if (j != k) {
                     sum2 += m_molefracs[j] / m_bdiff(j,k);
@@ -297,4 +297,67 @@ void GasTransport::getMixDiffCoeffs(doublereal* const d)
     }
 }
 
+void GasTransport::getMixDiffCoeffsMole(doublereal* const d)
+{
+    update_T();
+    update_C();
+
+    // update the binary diffusion coefficients if necessary
+    if (!m_bindiff_ok) {
+        updateDiff_T();
+    }
+
+    doublereal p = m_thermo->pressure();
+    if (m_nsp == 1) {
+        d[0] = m_bdiff(0,0) / p;
+    } else {
+        for (size_t k = 0; k < m_nsp; k++) {
+            double sum2 = 0.0;
+            for (size_t j = 0; j < m_nsp; j++) {
+                if (j != k) {
+                    sum2 += m_molefracs[j] / m_bdiff(j,k);
+                }
+            }
+            if (sum2 <= 0.0) {
+                d[k] = m_bdiff(k,k) / p;
+            } else {
+                d[k] = (1 - m_molefracs[k]) / (p * sum2);
+            }
+        }
+    }
+}
+
+void GasTransport::getMixDiffCoeffsMass(doublereal* const d)
+{
+    update_T();
+    update_C();
+
+    // update the binary diffusion coefficients if necessary
+    if (!m_bindiff_ok) {
+        updateDiff_T();
+    }
+
+    doublereal mmw = m_thermo->meanMolecularWeight();
+    doublereal p = m_thermo->pressure();
+
+    if (m_nsp == 1) {
+        d[0] = m_bdiff(0,0) / p;
+    } else {
+        for (size_t k=0; k<m_nsp; k++) {
+            double sum1 = 0.0;
+            double sum2 = 0.0;
+            for (size_t i=0; i<m_nsp; i++) {
+                if (i==k) {
+                    continue;
+                }
+                sum1 += m_molefracs[i] / m_bdiff(k,i);
+                sum2 += m_molefracs[i] * m_mw[i] / m_bdiff(k,i);
+            }
+            sum1 *= p;
+            sum2 *= p * m_molefracs[k] / (mmw - m_mw[k]*m_molefracs[k]);
+            d[k] = 1.0 / (sum1 +  sum2);
+        }
+    }
+}
+
 }