[1D] Add option for mass fraction gradient fluxes

include/cantera/oneD/StFlow.h

@@ -106,6 +106,24 @@ class StFlow : public Domain1D
         return m_do_soret;
     }
 
+    //! Compute species diffusive fluxes with respect to
+    //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass)
+    //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default)
+    //! when using the mixture-averaged transport model.
+    //! @param fluxGradientBasis set flux computation to mass or mole basis
+    //! @since New in %Cantera 3.1.
+    void setFluxGradientBasis(ThermoBasis fluxGradientBasis);
+
+    //! Compute species diffusive fluxes with respect to
+    //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass)
+    //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default)
+    //! when using the mixture-averaged transport model.
+    //! @return the basis used for flux computation (mass or mole fraction gradients)
+    //! @since New in %Cantera 3.1.
+    ThermoBasis fluxGradientBasis() const {
+        return m_fluxGradientBasis;
+    }
+
     //! Set the pressure. Since the flow equations are for the limit of small
     //! Mach number, the pressure is very nearly constant throughout the flow.
     void setPressure(double p) {
@@ -639,6 +657,7 @@ class StFlow : public Domain1D
     // flags
     vector<bool> m_do_energy;
     bool m_do_soret = false;
+    ThermoBasis m_fluxGradientBasis = ThermoBasis::molar;
     vector<bool> m_do_species;
     bool m_do_multicomponent = false;
 

interfaces/cython/cantera/_onedim.pxd

@@ -88,6 +88,8 @@ cdef extern from "cantera/oneD/StFlow.h":
         cbool doEnergy(size_t)
         void enableSoret(cbool) except +translate_exception
         cbool withSoret()
+        void setFluxGradientBasis(ThermoBasis) except +translate_exception
+        ThermoBasis fluxGradientBasis()
         void setFreeFlow()
         void setAxisymmetricFlow()
 

interfaces/cython/cantera/_onedim.pyx

@@ -515,6 +515,27 @@ cdef class _FlowBase(Domain1D):
         def __set__(self, enable):
             self.flow.enableSoret(<cbool>enable)
 
+    property flux_gradient_basis:
+        """
+        Get/Set whether or not species diffusive fluxes are computed with
+        respect to their mass fraction gradients ('mass')
+        or mole fraction gradients ('molar', default) when
+        using the mixture-averaged transport model.
+        """
+        def __get__(self):
+            if self.flow.fluxGradientBasis() == ThermoBasis.molar:
+                return 'molar'
+            else:
+                return 'mass'
+        def __set__(self, basis):
+            if basis == 'molar':
+                self.flow.setFluxGradientBasis(ThermoBasis.molar)
+            elif basis == 'mass':
+                self.flow.setFluxGradientBasis(ThermoBasis.mass)
+            else:
+                raise ValueError("Valid choices are 'mass' or 'molar'."
+                                 " Got {!r}.".format(basis))
+
     property energy_enabled:
         """ Determines whether or not to solve the energy equation."""
         def __get__(self):

interfaces/cython/cantera/onedim.py

@@ -215,6 +215,20 @@ def soret_enabled(self):
     def soret_enabled(self, enable):
         self.flame.soret_enabled = enable
 
+    @property
+    def flux_gradient_basis(self):
+        """
+        Get/Set whether or not species diffusive fluxes are computed with
+        respect to their mass fraction gradients ('mass')
+        or mole fraction gradients ('molar', default) when
+        using the mixture-averaged transport model.
+        """
+        return self.flame.flux_gradient_basis
+
+    @flux_gradient_basis.setter
+    def flux_gradient_basis(self, basis):
+        self.flame.flux_gradient_basis = basis
+
     @property
     def radiation_enabled(self):
         """

src/oneD/StFlow.cpp

@@ -217,6 +217,16 @@ string StFlow::transportModel() const {
     return m_trans->transportModel();
 }
 
+void StFlow::setFluxGradientBasis(ThermoBasis fluxGradientBasis) {
+    m_fluxGradientBasis = fluxGradientBasis;
+    if (transportModel() != "mixture-averaged-CK" &&
+        transportModel() != "mixture-averaged") {
+        warn_user("StFlow::setFluxGradientBasis",
+                  "Setting fluxGradientBasis only affects "
+                  "the mixture-averaged diffusion model.");
+    }
+}
+
 void StFlow::_getInitialSoln(double* x)
 {
     for (size_t j = 0; j < m_points; j++) {
@@ -629,11 +639,24 @@ void StFlow::updateTransport(double* x, size_t j0, size_t j1)
         for (size_t j = j0; j < j1; j++) {
             setGasAtMidpoint(x,j);
             m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
-            m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
+
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
+            } else {
+                m_trans->getMixDiffCoeffsMass(&m_diff[j*m_nsp]);
+            }
+
             double rho = m_thermo->density();
-            double wtm = m_thermo->meanMolecularWeight();
-            for (size_t k=0; k < m_nsp; k++) {
-                m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm;
+
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                double wtm = m_thermo->meanMolecularWeight();
+                for (size_t k=0; k < m_nsp; k++) {
+                    m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm;
+                }
+            } else {
+                for (size_t k=0; k < m_nsp; k++) {
+                    m_diff[k+j*m_nsp] *= rho;
+                }
             }
             m_tcon[j] = m_trans->thermalConductivity();
         }
@@ -674,10 +697,16 @@ void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
         for (size_t j = j0; j < j1; j++) {
             double sum = 0.0;
             double dz = z(j+1) - z(j);
-            for (size_t k = 0; k < m_nsp; k++) {
-                m_flux(k,j) = m_diff[k+m_nsp*j];
-                m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
-                sum -= m_flux(k,j);
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                for (size_t k = 0; k < m_nsp; k++) {
+                    m_flux(k,j) = m_diff[k+m_nsp*j] * (X(x,k,j) - X(x,k,j+1))/dz;
+                    sum -= m_flux(k,j);
+                }
+            } else {
+                for (size_t k = 0; k < m_nsp; k++) {
+                    m_flux(k,j) = m_diff[k+m_nsp*j] * (Y(x,k,j) - Y(x,k,j+1))/dz;
+                    sum -= m_flux(k,j);
+                }
             }
             // correction flux to insure that \sum_k Y_k V_k = 0.
             for (size_t k = 0; k < m_nsp; k++) {
@@ -780,6 +809,8 @@ AnyMap StFlow::getMeta() const
 
     state["Soret-enabled"] = m_do_soret;
 
+    state["flux-gradient-basis"] = static_cast<long int>(m_fluxGradientBasis);
+
     set<bool> species_flags(m_do_species.begin(), m_do_species.end());
     if (species_flags.size() == 1) {
         state["species-enabled"] = m_do_species[0];
@@ -886,6 +917,11 @@ void StFlow::setMeta(const AnyMap& state)
         m_do_soret = state["Soret-enabled"].asBool();
     }
 
+    if (state.hasKey("flux-gradient-basis")) {
+        m_fluxGradientBasis = static_cast<ThermoBasis>(
+                state["flux-gradient-basis"].asInt());
+    }
+
     if (state.hasKey("species-enabled")) {
         const AnyValue& se = state["species-enabled"];
         if (se.isScalar()) {

test/python/test_onedim.py

@@ -480,6 +480,22 @@ def test_unity_lewis(self):
         # towards zero as grid is refined)
         self.assertLess(dh_unity_lewis, 0.1 * dh_mix)
 
+    def test_flux_gradient_mass(self):
+        self.create_sim(ct.one_atm, 300, 'H2:1.1, O2:1, AR:5.3')
+        self.sim.transport_model = 'mixture-averaged'
+        self.sim.set_refine_criteria(ratio=3.0, slope=0.08, curve=0.12)
+        assert self.sim.flux_gradient_basis == 'molar'
+        self.sim.solve(loglevel=0, auto=True)
+        sl_mole = self.sim.velocity[0]
+        self.sim.flux_gradient_basis = 'mass'
+        assert self.sim.flux_gradient_basis == 'mass'
+        self.sim.solve(loglevel=0)
+        sl_mass = self.sim.velocity[0]
+        # flame speeds should not be exactly equal
+        assert sl_mole != sl_mass
+        # but they should be close
+        assert sl_mole == approx(sl_mass, rel=0.1)
+
     def test_soret_with_mix(self):
         # Test that enabling Soret diffusion without
         # multicomponent transport results in an error