[1D] Remove charge neutrality solver option from IonFlow

lightsquared · Jun 15, 2018 · eeb27d8 · eeb27d8
1 parent 3a0f46e
commit eeb27d8
Show file tree

Hide file tree

Showing 7 changed files with 10 additions and 201 deletions.
diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h
@@ -18,16 +18,7 @@ namespace Cantera
  * diffusion rate of electron without internal electric forces (ambi-
  * polar diffusion effect).
  *
- * The second stage uses charge neutrality model, which assume zero charge
- * flux throughout the domain, to calculate drift flux. The drift flux is
- * added to the total flux of ions.
- * Reference:
- * Prager, J., U. Riedel, and J. Warnatz.
- * "Modeling ion chemistry and charged species diffusion in lean
- * methane–oxygen flames."
- * Proceedings of the Combustion Institute 31.1 (2007): 1129-1137.
- *
- * The third stage evaluates drift flux from electric field calculated from
+ * The second stage evaluates drift flux from electric field calculated from
  * Poisson's equation, which is solved together with other equations. Poisson's
  * equation is coupled because the total charge densities depends on the species'
  * concentration.
@@ -56,13 +47,6 @@ class IonFlow : public FreeFlame
     bool doPoisson(size_t j) {
         return m_do_poisson[j];
     }
-    //! set to solve velocity on a point
-    void solveVelocity(size_t j=npos);
-    //! set to fix velocity on a point
-    void fixVelocity(size_t j=npos);
-    bool doVelocity(size_t j) {
-        return m_do_velocity[j];
-    }
 
     /**
      * Sometimes it is desired to carry out the simulation using a specified
@@ -89,14 +73,10 @@ class IonFlow : public FreeFlame
     virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1);
     //! Solving phase one: the fluxes of charged species are turned off
     virtual void frozenIonMethod(const double* x, size_t j0, size_t j1);
-    //! Solving phase two: the Prager's ambipolar-diffusion model is used
-    virtual void chargeNeutralityModel(const double* x, size_t j0, size_t j1);
     //! Solving phase three: the Poisson's equation is added coupled by the electrical drift
     virtual void poissonEqnMethod(const double* x, size_t j0, size_t j1);
     //! flag for solving poisson's equation or not
     std::vector<bool> m_do_poisson;
-    //! flag for solving the velocity or not
-    std::vector<bool> m_do_velocity;
 
     //! flag for importing transport of electron
     bool m_import_electron_transport;
@@ -133,9 +113,6 @@ class IonFlow : public FreeFlame
     //! fixed electric potential value
     vector_fp m_fixedElecPoten;
 
-    //! fixed velocity value
-    vector_fp m_fixedVelocity;
-
     //! fixed electron transport values
     vector_fp m_ztfix;
     vector_fp m_diff_e_fix;
@@ -146,11 +123,6 @@ class IonFlow : public FreeFlame
         return m_fixedElecPoten[j];
     }
 
-    //! The fixed velocity value at point j
-    double u_fixed(size_t j) const {
-        return m_fixedVelocity[j];
-    }
-
     //! electric potential
     double phi(const double* x, size_t j) const {
         return x[index(c_offset_P, j)];

diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
@@ -703,9 +703,6 @@ cdef extern from "cantera/oneD/IonFlow.h":
         void solvePoissonEqn()
         void fixElectricPotential()
         cbool doPoisson(size_t)
-        void solveVelocity()
-        void fixVelocity()
-        cbool doVelocity(size_t)
 
 
 cdef extern from "cantera/oneD/Sim1D.h":

diff --git a/interfaces/cython/cantera/examples/onedim/ion_flame.py b/interfaces/cython/cantera/examples/onedim/ion_flame.py
@@ -19,19 +19,15 @@
 
 # Set up flame object
 f = ct.IonFlame(gas, width=width)
-f.set_refine_criteria(ratio=3, slope=0.06, curve=0.12)
+f.set_refine_criteria(ratio=3, slope=0.05, curve=0.1)
 f.show_solution()
 
 # stage one
 f.solve(loglevel=loglevel, auto=True)
 
 # stage two
-f.solve(loglevel=loglevel, stage=2, enable_energy=False)
 f.solve(loglevel=loglevel, stage=2, enable_energy=True)
 
-# stage three
-f.solve(loglevel=loglevel, stage=3, enable_energy=True)
-
 f.save('CH4_adiabatic.xml', 'mix', 'solution with mixture-averaged transport')
 f.show_solution()
 print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))

diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
@@ -570,20 +570,10 @@ def __init__(self, gas, grid=None, width=None):
         super(IonFlame, self).__init__(gas, grid, width)
 
     def solve(self, loglevel=1, refine_grid=True, auto=False, stage=1, enable_energy=True):
-        if enable_energy == True:
-            self.energy_enabled = True
-            self.velocity_enabled = True
-        else:
-            self.energy_enabled = False
-            self.velocity_enabled = False
+        self.flame.set_solvingStage(stage)
         if stage == 1:
-            self.flame.set_solvingStage(stage)
             super(IonFlame, self).solve(loglevel, refine_grid, auto)
         if stage == 2:
-            self.flame.set_solvingStage(stage)
-            super(IonFlame, self).solve(loglevel, refine_grid, auto)
-        if stage == 3:
-            self.flame.set_solvingStage(stage)
             self.poisson_enabled = True
             super(IonFlame, self).solve(loglevel, refine_grid, auto)
 
@@ -626,15 +616,6 @@ def poisson_enabled(self):
     def poisson_enabled(self, enable):
         self.flame.poisson_enabled = enable
 
-    @property
-    def velocity_enabled(self):
-        """ Get/Set whether or not to solve the velocity."""
-        return self.flame.velocity_enabled
-
-    @velocity_enabled.setter
-    def velocity_enabled(self, enable):
-        self.flame.velocity_enabled = enable
-
     @property
     def phi(self):
         """

diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx
@@ -506,16 +506,6 @@ cdef class IonFlow(_FlowBase):
             else:
                 (<CxxIonFlow*>self.flow).fixElectricPotential()
 
-    property velocity_enabled:
-        """ Determines whether or not to solve the velocity."""
-        def __get__(self):
-            return (<CxxIonFlow*>self.flow).doVelocity(0)
-        def __set__(self, enable):
-            if enable:
-                (<CxxIonFlow*>self.flow).solveVelocity()
-            else:
-                (<CxxIonFlow*>self.flow).fixVelocity()
-
 
 cdef class AxisymmetricStagnationFlow(_FlowBase):
     """

diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py
@@ -946,14 +946,8 @@ def test_ion_profile(self):
         # stage one
         self.sim.solve(loglevel=0, auto=True)
 
-        # stage two
-        self.sim.solve(loglevel=0, stage=2, enable_energy=False)
-
-        # stage two
+        #stage two
         self.sim.solve(loglevel=0, stage=2, enable_energy=True)
 
-        #stage three
-        self.sim.solve(loglevel=0, stage=3, enable_energy=True)
-
         # Regression test
-        self.assertNear(min(self.sim.E) / max(self.sim.E), -5.0765, 1e-3)
+        self.assertNear(max(self.sim.E), 113.5274, 1e-3)
diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp
@@ -48,17 +48,14 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     m_refiner->setActive(c_offset_P, false);
     m_mobility.resize(m_nsp*m_points);
     m_do_poisson.resize(m_points,false);
-    m_do_velocity.resize(m_points,true);
 }
 
 void IonFlow::resize(size_t components, size_t points){
     StFlow::resize(components, points);
     m_mobility.resize(m_nsp*m_points);
     m_do_species.resize(m_nsp,true);
     m_do_poisson.resize(m_points,false);
-    m_do_velocity.resize(m_points,true);
     m_fixedElecPoten.resize(m_points,0.0);
-    m_fixedVelocity.resize(m_points);
     m_elecMobility.resize(m_points);
     m_elecDiffCoeff.resize(m_points);
 }
@@ -87,9 +84,6 @@ void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
         frozenIonMethod(x,j0,j1);
     }
     if (m_stage == 2) {
-        chargeNeutralityModel(x,j0,j1);
-    }
-    if (m_stage == 3) {
         poissonEqnMethod(x,j0,j1);
     }
 }
@@ -121,55 +115,6 @@ void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1)
     }
 }
 
-void IonFlow::chargeNeutralityModel(const double* x, size_t j0, size_t j1)
-{
-    for (size_t j = j0; j < j1; j++) {
-        double wtm = m_wtm[j];
-        double rho = density(j);
-        double dz = z(j+1) - z(j);
-
-        // mixture-average diffusion
-        for (size_t k = 0; k < m_nsp; k++) {
-            m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
-            m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
-        }
-
-        // ambipolar diffusion
-        double sum_chargeFlux = 0.0;
-        double sum = 0.0;
-        for (size_t k : m_kCharge) {
-            double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j));
-            int q_k = m_speciesCharge[k];
-            sum_chargeFlux += m_speciesCharge[k] / m_wt[k] * m_flux(k,j);
-            // The mobility is used because it is more general than
-            // using diffusion coefficient and Einstein relation
-            sum += m_mobility[k+m_nsp*j] * Xav * q_k * q_k;
-        }
-        for (size_t k : m_kCharge) {
-            double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j));
-            double drift;
-            int q_k = m_speciesCharge[k];
-            drift = q_k * q_k * m_mobility[k+m_nsp*j] * Xav / sum;
-            drift *= -sum_chargeFlux * m_wt[k] / q_k;
-            m_flux(k,j) += drift;
-        }
-
-        // correction flux
-        double sum_flux = 0.0;
-        for (size_t k = 0; k < m_nsp; k++) {
-            sum_flux -= m_flux(k,j); // total net flux
-        }
-        double sum_ion = 0.0;
-        for (size_t k : m_kCharge) {
-            sum_ion += Y(x,k,j);
-        }
-        // The portion of correction for ions is taken off
-        for (size_t k : m_kNeutral) {
-            m_flux(k,j) += Y(x,k,j) / (1-sum_ion) * sum_flux;
-        }
-    }
-}
-
 void IonFlow::poissonEqnMethod(const double* x, size_t j0, size_t j1)
 {
     for (size_t j = j0; j < j1; j++) {
@@ -212,14 +157,13 @@ void IonFlow::poissonEqnMethod(const double* x, size_t j0, size_t j1)
 
 void IonFlow::setSolvingStage(const size_t stage)
 {
-    if (stage == 1 || stage == 2 || stage == 3) {
+    if (stage == 1 || stage == 2) {
         m_stage = stage;
     } else {
         throw CanteraError("IonFlow::updateDiffFluxes",
-                    "solution phase must be set to:"
-                    "1: frozenIonMethod"
-                    "2: chargeNeutralityModel"
-                    "3: poissonEqnMethod");
+                    "solution stage must be set to: "
+                    "1) frozenIonMethod, "
+                    "2) poissonEqnMethod");
     }
 }
 
@@ -234,7 +178,7 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag,
                            double rdt, size_t jmin, size_t jmax)
 {
     StFlow::evalResidual(x, rsd, diag, rdt, jmin, jmax);
-    if (m_stage != 3) {
+    if (m_stage != 2) {
         return;
     }
 
@@ -259,13 +203,6 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag,
             }
             rsd[index(c_offset_P, j)] = dEdz(x,j) - chargeDensity / epsilon_0;
             diag[index(c_offset_P, j)] = 0;
-
-            // This method is used when you disable energy equation
-            // but still maintain the velocity profile
-            if (!m_do_velocity[j]) {
-                rsd[index(c_offset_U, j)] = u(x,j) - u_fixed(j);
-                diag[index(c_offset_U, j)] = 0;
-            }
         }
     }
 }
@@ -320,54 +257,6 @@ void IonFlow::fixElectricPotential(size_t j)
     }
 }
 
-void IonFlow::solveVelocity(size_t j)
-{
-    bool changed = false;
-    if (j == npos) {
-        for (size_t i = 0; i < m_points; i++) {
-            if (!m_do_velocity[i]) {
-                changed = true;
-            }
-            m_do_velocity[i] = true;
-        }
-    } else {
-        if (!m_do_velocity[j]) {
-            changed = true;
-        }
-        m_do_velocity[j] = true;
-    }
-    m_refiner->setActive(c_offset_U, true);
-    m_refiner->setActive(c_offset_V, true);
-    m_refiner->setActive(c_offset_T, true);
-    if (changed) {
-        needJacUpdate();
-    }
-}
-
-void IonFlow::fixVelocity(size_t j)
-{
-    bool changed = false;
-    if (j == npos) {
-        for (size_t i = 0; i < m_points; i++) {
-            if (m_do_velocity[i]) {
-                changed = true;
-            }
-            m_do_velocity[i] = false;
-        }
-    } else {
-        if (m_do_velocity[j]) {
-            changed = true;
-        }
-        m_do_velocity[j] = false;
-    }
-    m_refiner->setActive(c_offset_U, false);
-    m_refiner->setActive(c_offset_V, false);
-    m_refiner->setActive(c_offset_T, false);
-    if (changed) {
-        needJacUpdate();
-    }
-}
-
 void IonFlow::setElectronTransport(vector_fp& zfixed, vector_fp& diff_e_fixed,
                                    vector_fp& mobi_e_fixed)
 {
@@ -390,16 +279,6 @@ void IonFlow::_finalize(const double* x)
     if (p) {
         solvePoissonEqn();
     }
-    // save the velocity profile if the velocity is disabled
-    bool v = m_do_velocity[0];
-    for (size_t j = 0; j < m_points; j++) {
-        if (!v) {
-            m_fixedVelocity[j] = u(x,j);
-        }
-    }
-    if (v) {
-        solveVelocity();
-    }
 }
 
 }