turbulent diffusivity added to increase thermal diffusivity

include/cantera/transport/ConstantTransport.h

@@ -0,0 +1,52 @@
+/**
+ *  @file ConstantTransport.h
+ *    Headers for the ConstantTransport object, which models transport
+ *    properties in ideal gas solutions using the unity Lewis number
+ *    approximation
+ *    (see \ref tranprops and \link Cantera::ConstantTransport ConstantTransport \endlink) .
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_CONSTANTTRAN_H
+#define CT_CONSTANTTRAN_H
+
+#include "MixTransport.h"
+
+namespace Cantera
+{
+//! Class ConstantTransport implements the unity Lewis number approximation
+//! for the mixture-averaged species diffusion coefficients. Mixture-averaged
+//! transport properties for viscosity and thermal conductivity are inherited
+//! from the MixTransport class.
+//! @ingroup tranprops
+class ConstantTransport : public MixTransport
+{
+public:
+    //! Default constructor.
+    ConstantTransport();
+
+    virtual std::string transportType() const {
+        return "Constant";
+    }
+
+    //! Update the internal parameters whenever the temperature has changed
+    /*!
+     * This is called whenever a transport property is requested if the
+     * temperature has changed since the last call to update_T().
+     */
+    virtual void update_T();
+
+    virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
+
+protected:
+    //! Reference temperature at which binary diffusivities are calculated
+    /*!
+     *  Units = K
+     */
+    doublereal m_reftemp;
+
+};
+}
+#endif

include/cantera/transport/TurbulentTransport.h

@@ -0,0 +1,121 @@
+/**
+ *  @file TurbulentTransport.h
+ *    Headers for the TurbulentTransport object, which models transport
+ *    properties in ideal gas solutions using the unity Lewis number
+ *    approximation
+ *    (see \ref tranprops and \link Cantera::TurbulentTransport TurbulentTransport \endlink) .
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_TURBULENTTRAN_H
+#define CT_TURBULENTTRAN_H
+
+#include "MixTransport.h"
+
+namespace Cantera
+{
+//! Class TurbulentTransport implements the unity Lewis number approximation
+//! for the mixture-averaged species diffusion coefficients. Mixture-averaged
+//! transport properties for viscosity and thermal conductivity are inherited
+//! from the MixTransport class.
+//! @ingroup tranprops
+class TurbulentTransport : public MixTransport
+{
+public:
+    TurbulentTransport();
+
+    virtual std::string transportType() const {
+        return "Turbulent";
+    }
+
+    virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
+
+    //! Returns the mixture thermal conductivity (W/m /K)
+    /*!
+     * The thermal conductivity is computed from the following mixture rule:
+     * \f[
+     *     \lambda = 0.5 \left( \sum_k X_k \lambda_k  + \frac{1}{\sum_k X_k/\lambda_k} \right)
+     * \f]
+     *
+     * It's used to compute the flux of energy due to a thermal gradient
+     *
+     * \f[
+     *     j_T =  - \lambda  \nabla T
+     * \f]
+     *
+     * The flux of energy has units of energy (kg m2 /s2) per second per area.
+     *
+     * The units of lambda are W / m K which is equivalent to kg m / s^3 K.
+     *
+     * @returns the mixture thermal conductivity, with units of W/m/K
+     */
+    virtual doublereal thermalConductivity();
+
+    //! Returns the unity Lewis number approximation based diffusion
+    //! coefficients [m^2/s].
+    /*!
+     * Returns the unity Lewis number approximation based diffusion coefficients
+     * for a gas, appropriate for calculating the mass averaged diffusive flux
+     * with respect to the mass averaged velocity using gradients of the mole
+     * fraction.
+     *
+     * \f[
+     *     D^\prime_{km} = \frac{\lambda}{\rho c_p}
+     * \f]
+     *
+     * In order to obtain the expected behavior from a unity Lewis number model,
+     * this formulation requires that the correction velocity be computed as
+     *
+     * \f[
+     *     V_c = \sum \frac{W_k}{\overline{W}} D^\prime_{km} \nabla X_k
+     * \f]
+     *
+     * @param[out] d  Vector of diffusion coefficients for each species (m^2/s).
+     * length m_nsp.
+     */
+    virtual void getMixDiffCoeffs(double* const d) {
+        GasTransport::getMixDiffCoeffs(d);
+
+        for (size_t k = 0; k < m_nsp; k++) {
+            d[k] += m_turbmodifier;
+        }
+    }
+
+    //! Not implemented for unity Lewis number approximation
+    virtual void getMixDiffCoeffsMole(double* const d){
+        throw NotImplementedError("TurbulentTransport::getMixDiffCoeffsMole");
+    }
+
+    //! Returns the unity Lewis number approximation based diffusion
+    //! coefficients [m^2/s].
+    /*!
+     * These are the coefficients for calculating the diffusive mass fluxes
+     * from the species mass fraction gradients, computed as
+     *
+     * \f[
+     *     D_{km} = \frac{\lambda}{\rho c_p}
+     * \f]
+     *
+     * @param[out] d  Vector of diffusion coefficients for each species (m^2/s).
+     * length m_nsp.
+     */
+    virtual void getMixDiffCoeffsMass(double* const d){
+        GasTransport::getMixDiffCoeffsMass(d);
+
+        for (size_t k = 0; k < m_nsp; k++) {
+            d[k] += m_turbmodifier;
+        }
+    }
+
+protected:
+
+    //! Internal storage for weights to calculate mean diffusivity
+    doublereal m_turbmodifier;
+    doublereal m_profwidth;
+    doublereal m_rdlayerstart;
+
+};
+}
+#endif

include/cantera/transport/WmeanTransport.h

@@ -0,0 +1,106 @@
+/**
+ *  @file WmeanTransport.h
+ *    Headers for the WmeanTransport object, which models transport
+ *    properties in ideal gas solutions using the unity Lewis number
+ *    approximation
+ *    (see \ref tranprops and \link Cantera::WmeanTransport WmeanTransport \endlink) .
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_WMEANTRAN_H
+#define CT_WMEANTRAN_H
+
+#include "MixTransport.h"
+
+namespace Cantera
+{
+//! Class WmeanTransport implements the unity Lewis number approximation
+//! for the mixture-averaged species diffusion coefficients. Mixture-averaged
+//! transport properties for viscosity and thermal conductivity are inherited
+//! from the MixTransport class.
+//! @ingroup tranprops
+class WmeanTransport : public MixTransport
+{
+public:
+    WmeanTransport();
+
+    virtual std::string transportType() const {
+        return "Wmean";
+    }
+
+    virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
+
+    //! Returns the unity Lewis number approximation based diffusion
+    //! coefficients [m^2/s].
+    /*!
+     * Returns the unity Lewis number approximation based diffusion coefficients
+     * for a gas, appropriate for calculating the mass averaged diffusive flux
+     * with respect to the mass averaged velocity using gradients of the mole
+     * fraction.
+     *
+     * \f[
+     *     D^\prime_{km} = \frac{\lambda}{\rho c_p}
+     * \f]
+     *
+     * In order to obtain the expected behavior from a unity Lewis number model,
+     * this formulation requires that the correction velocity be computed as
+     *
+     * \f[
+     *     V_c = \sum \frac{W_k}{\overline{W}} D^\prime_{km} \nabla X_k
+     * \f]
+     *
+     * @param[out] d  Vector of diffusion coefficients for each species (m^2/s).
+     * length m_nsp.
+     */
+    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];
+        }
+        for (size_t k = 0; k < m_nsp; k++) {
+            d[k] = Dm;
+        }
+    }
+
+    //! Not implemented for unity Lewis number approximation
+    virtual void getMixDiffCoeffsMole(double* const d){
+        throw NotImplementedError("WmeanTransport::getMixDiffCoeffsMole");
+    }
+
+    //! Returns the unity Lewis number approximation based diffusion
+    //! coefficients [m^2/s].
+    /*!
+     * These are the coefficients for calculating the diffusive mass fluxes
+     * from the species mass fraction gradients, computed as
+     *
+     * \f[
+     *     D_{km} = \frac{\lambda}{\rho c_p}
+     * \f]
+     *
+     * @param[out] d  Vector of diffusion coefficients for each species (m^2/s).
+     * length m_nsp.
+     */
+    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];
+        }
+        for (size_t k = 0; k < m_nsp; k++) {
+            d[k] = Dm;
+        }
+    }
+
+protected:
+
+    //! Internal storage for weights to calculate mean diffusivity
+    vector_fp m_weight;
+
+};
+}
+#endif

src/transport/ConstantTransport.cpp

@@ -0,0 +1,162 @@
+/**
+ *  @file ConstantTransport.cpp
+ *  Mixture-averaged transport properties for ideal gas mixtures.
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/transport/ConstantTransport.h"
+#include "cantera/base/stringUtils.h"
+
+using namespace std;
+
+namespace Cantera
+{
+ConstantTransport::ConstantTransport() :
+    MixTransport(),
+    m_reftemp(300.0)
+{
+}
+
+void ConstantTransport::init(ThermoPhase* thermo, int mode, int log_level)
+{
+    MixTransport::init(thermo, mode, log_level);
+    // m_reftemp = thermo->temperature();
+    // cout << "check reference temperature" << m_reftemp << endl;
+}
+
+    /*
+void ConstantTransport::init(ThermoPhase* thermo, int mode, int log_level)
+{
+    MixTransport::init(thermo, mode, log_level);
+}
+
+void ConstantTransport::getMobilities(doublereal* const mobil)
+{
+    getMixDiffCoeffs(m_spwork.data());
+    doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
+    for (size_t k = 0; k < m_nsp; k++) {
+        mobil[k] = c1 * m_spwork[k];
+    }
+}
+
+doublereal ConstantTransport::thermalConductivity()
+{
+    update_T();
+    update_C();
+    if (!m_spcond_ok) {
+        updateCond_T();
+    }
+    if (!m_condmix_ok) {
+        doublereal sum1 = 0.0, sum2 = 0.0;
+        for (size_t k = 0; k < m_nsp; k++) {
+            sum1 += m_molefracs[k] * m_cond[k];
+            sum2 += m_molefracs[k] / m_cond[k];
+        }
+        m_lambda = 0.5*(sum1 + 1.0/sum2);
+        m_condmix_ok = true;
+    }
+    return m_lambda;
+}
+
+void ConstantTransport::getThermalDiffCoeffs(doublereal* const dt)
+{
+    for (size_t k = 0; k < m_nsp; k++) {
+        dt[k] = 0.0;
+    }
+}
+
+void ConstantTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
+                                    size_t ldx, const doublereal* const grad_X,
+                                    size_t ldf, doublereal* const fluxes)
+{
+    update_T();
+    update_C();
+    getMixDiffCoeffs(m_spwork.data());
+    const vector_fp& mw = m_thermo->molecularWeights();
+    const doublereal* y = m_thermo->massFractions();
+    doublereal rhon = m_thermo->molarDensity();
+    vector_fp sum(ndim,0.0);
+    for (size_t n = 0; n < ndim; n++) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k];
+            sum[n] += fluxes[n*ldf + k];
+        }
+    }
+    // add correction flux to enforce sum to zero
+    for (size_t n = 0; n < ndim; n++) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            fluxes[n*ldf + k] -= y[k]*sum[n];
+        }
+    }
+}
+
+void ConstantTransport::update_T()
+{
+    doublereal t = m_thermo->temperature();
+    if (t == m_temp && m_nsp == m_thermo->nSpecies()) {
+        return;
+    }
+    if (t < 0.0) {
+        throw CanteraError("ConstantTransport::update_T",
+                           "negative temperature {}", t);
+    }
+    GasTransport::update_T();
+    // temperature has changed, so polynomial fits will need to be redone.
+    m_spcond_ok = false;
+    m_bindiff_ok = false;
+    m_condmix_ok = false;
+}
+
+void ConstantTransport::update_C()
+{
+    // signal that concentration-dependent quantities will need to be recomputed
+    // before use, and update the local mole fractions.
+    m_visc_ok = false;
+    m_condmix_ok = false;
+    m_thermo->getMoleFractions(m_molefracs.data());
+
+    // add an offset to avoid a pure species condition
+    for (size_t k = 0; k < m_nsp; k++) {
+        m_molefracs[k] = std::max(Tiny, m_molefracs[k]);
+    }
+}
+
+void ConstantTransport::updateCond_T()
+{
+    if (m_mode == CK_Mode) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k]));
+        }
+    } else {
+        for (size_t k = 0; k < m_nsp; k++) {
+            m_cond[k] = m_sqrt_t * dot5(m_polytempvec, m_condcoeffs[k]);
+        }
+    }
+    m_spcond_ok = true;
+    m_condmix_ok = false;
+}
+
+    */
+
+void ConstantTransport::update_T()
+{
+    doublereal t = m_thermo->temperature();
+    if (m_reftemp == m_temp && m_nsp == m_thermo->nSpecies()) {
+        return;
+    }
+    if (t < 0.0) {
+        throw CanteraError("ConstantTransport::update_T",
+                           "negative temperature {}", t);
+    }
+    m_thermo->setTemperature(m_reftemp);
+    GasTransport::update_T();
+    m_thermo->setTemperature(t);
+    // temperature has changed, so polynomial fits will need to be redone.
+    m_spcond_ok = false;
+    m_bindiff_ok = false;
+    m_condmix_ok = false;
+}
+
+}

src/transport/TransportFactory.cpp

@@ -7,6 +7,9 @@
 #include "cantera/transport/MultiTransport.h"
 #include "cantera/transport/MixTransport.h"
 #include "cantera/transport/UnityLewisTransport.h"
+#include "cantera/transport/ConstantTransport.h"
+#include "cantera/transport/WmeanTransport.h"
+#include "cantera/transport/TurbulentTransport.h"
 #include "cantera/transport/IonGasTransport.h"
 #include "cantera/transport/WaterTransport.h"
 #include "cantera/transport/DustyGasTransport.h"
@@ -45,6 +48,12 @@ TransportFactory::TransportFactory()
 {
     reg("", []() { return new Transport(); });
     m_synonyms["None"] = "";
+    reg("Constant", []() { return new ConstantTransport(); });
+    m_synonyms["constant"] = "Constant";
+    reg("Wmean", []() { return new WmeanTransport(); });
+    m_synonyms["wmean"] = "Wmean";
+    reg("Turbulent", []() { return new TurbulentTransport(); });
+    m_synonyms["turbulent"] = "Turbulent";
     reg("UnityLewis", []() { return new UnityLewisTransport(); });
     m_synonyms["unity-Lewis-number"] = "UnityLewis";
     reg("Mix", []() { return new MixTransport(); });

src/transport/TurbulentTransport.cpp

@@ -0,0 +1,49 @@
+/**
+ *  @file TurbulentTransport.cpp
+ *  Mixture-averaged transport properties for ideal gas mixtures.
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/transport/TurbulentTransport.h"
+#include "cantera/base/stringUtils.h"
+
+using namespace std;
+
+namespace Cantera
+{
+TurbulentTransport::TurbulentTransport() :
+    MixTransport(),
+    m_turbmodifier(1.0),
+    m_profwidth(1e-5),
+    m_rdlayerstart(300.)
+{
+}
+
+void TurbulentTransport::init(ThermoPhase* thermo, int mode, int log_level)
+{
+    MixTransport::init(thermo, mode, log_level);
+
+    ifstream myfile;
+    myfile.open ("turbulent-coefs.txt");
+
+    myfile >> m_turbmodifier;
+
+    myfile.close();
+
+    // m_reftemp = thermo->temperature();
+    // cout << "check reference temperature" << m_reftemp << endl;
+}
+
+doublereal TurbulentTransport::thermalConductivity()
+{
+    doublereal rho = m_thermo->density();
+    doublereal cp = m_thermo->cp_mass();
+
+    MixTransport::thermalConductivity();
+    m_lambda += m_turbmodifier * rho * cp;
+    m_condmix_ok = true;
+    return m_lambda;
+}
+}

src/transport/WmeanTransport.cpp

@@ -0,0 +1,177 @@
+/**
+ *  @file WmeanTransport.cpp
+ *  Mixture-averaged transport properties for ideal gas mixtures.
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/transport/WmeanTransport.h"
+#include "cantera/base/stringUtils.h"
+
+using namespace std;
+
+namespace Cantera
+{
+WmeanTransport::WmeanTransport() :
+    MixTransport(),
+    m_weight(m_nsp)
+{
+}
+
+void WmeanTransport::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];
+
+        cout << m_thermo->speciesName(k) << m_weight[k] << endl;
+    }
+    myfile.close();
+
+    // m_reftemp = thermo->temperature();
+    // cout << "check reference temperature" << m_reftemp << endl;
+}
+
+    /*
+void WmeanTransport::init(ThermoPhase* thermo, int mode, int log_level)
+{
+    MixTransport::init(thermo, mode, log_level);
+}
+
+void WmeanTransport::getMobilities(doublereal* const mobil)
+{
+    getMixDiffCoeffs(m_spwork.data());
+    doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
+    for (size_t k = 0; k < m_nsp; k++) {
+        mobil[k] = c1 * m_spwork[k];
+    }
+}
+
+doublereal WmeanTransport::thermalConductivity()
+{
+    update_T();
+    update_C();
+    if (!m_spcond_ok) {
+        updateCond_T();
+    }
+    if (!m_condmix_ok) {
+        doublereal sum1 = 0.0, sum2 = 0.0;
+        for (size_t k = 0; k < m_nsp; k++) {
+            sum1 += m_molefracs[k] * m_cond[k];
+            sum2 += m_molefracs[k] / m_cond[k];
+        }
+        m_lambda = 0.5*(sum1 + 1.0/sum2);
+        m_condmix_ok = true;
+    }
+    return m_lambda;
+}
+
+void WmeanTransport::getThermalDiffCoeffs(doublereal* const dt)
+{
+    for (size_t k = 0; k < m_nsp; k++) {
+        dt[k] = 0.0;
+    }
+}
+
+void WmeanTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
+                                    size_t ldx, const doublereal* const grad_X,
+                                    size_t ldf, doublereal* const fluxes)
+{
+    update_T();
+    update_C();
+    getMixDiffCoeffs(m_spwork.data());
+    const vector_fp& mw = m_thermo->molecularWeights();
+    const doublereal* y = m_thermo->massFractions();
+    doublereal rhon = m_thermo->molarDensity();
+    vector_fp sum(ndim,0.0);
+    for (size_t n = 0; n < ndim; n++) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k];
+            sum[n] += fluxes[n*ldf + k];
+        }
+    }
+    // add correction flux to enforce sum to zero
+    for (size_t n = 0; n < ndim; n++) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            fluxes[n*ldf + k] -= y[k]*sum[n];
+        }
+    }
+}
+
+void WmeanTransport::update_T()
+{
+    doublereal t = m_thermo->temperature();
+    if (t == m_temp && m_nsp == m_thermo->nSpecies()) {
+        return;
+    }
+    if (t < 0.0) {
+        throw CanteraError("WmeanTransport::update_T",
+                           "negative temperature {}", t);
+    }
+    GasTransport::update_T();
+    // temperature has changed, so polynomial fits will need to be redone.
+    m_spcond_ok = false;
+    m_bindiff_ok = false;
+    m_condmix_ok = false;
+}
+
+void WmeanTransport::update_C()
+{
+    // signal that concentration-dependent quantities will need to be recomputed
+    // before use, and update the local mole fractions.
+    m_visc_ok = false;
+    m_condmix_ok = false;
+    m_thermo->getMoleFractions(m_molefracs.data());
+
+    // add an offset to avoid a pure species condition
+    for (size_t k = 0; k < m_nsp; k++) {
+        m_molefracs[k] = std::max(Tiny, m_molefracs[k]);
+    }
+}
+
+void WmeanTransport::updateCond_T()
+{
+    if (m_mode == CK_Mode) {
+        for (size_t k = 0; k < m_nsp; k++) {
+            m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k]));
+        }
+    } else {
+        for (size_t k = 0; k < m_nsp; k++) {
+            m_cond[k] = m_sqrt_t * dot5(m_polytempvec, m_condcoeffs[k]);
+        }
+    }
+    m_spcond_ok = true;
+    m_condmix_ok = false;
+}
+
+void WmeanTransport::update_T()
+{
+    doublereal t = m_thermo->temperature();
+    if (m_reftemp == m_temp && m_nsp == m_thermo->nSpecies()) {
+        return;
+    }
+    if (t < 0.0) {
+        throw CanteraError("WmeanTransport::update_T",
+                           "negative temperature {}", t);
+    }
+    m_thermo->setTemperature(m_reftemp);
+    GasTransport::update_T();
+    m_thermo->setTemperature(t);
+    // temperature has changed, so polynomial fits will need to be redone.
+    m_spcond_ok = false;
+    m_bindiff_ok = false;
+    m_condmix_ok = false;
+}
+
+    */
+
+}