From eeb27d84a9c1a374da9bf7b4409370c62580aa84 Mon Sep 17 00:00:00 2001 From: bangshiuh Date: Fri, 18 May 2018 16:45:38 -0400 Subject: [PATCH] [1D] Remove charge neutrality solver option from IonFlow --- include/cantera/oneD/IonFlow.h | 30 +--- interfaces/cython/cantera/_cantera.pxd | 3 - .../cantera/examples/onedim/ion_flame.py | 6 +- interfaces/cython/cantera/onedim.py | 21 +-- interfaces/cython/cantera/onedim.pyx | 10 -- interfaces/cython/cantera/test/test_onedim.py | 10 +- src/oneD/IonFlow.cpp | 131 +----------------- 7 files changed, 10 insertions(+), 201 deletions(-) diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index 6143c2cee..39c457f02 100644 --- 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 @@ public: 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 @@ protected: 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 m_do_poisson; - //! flag for solving the velocity or not - std::vector m_do_velocity; //! flag for importing transport of electron bool m_import_electron_transport; @@ -133,9 +113,6 @@ protected: //! 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 @@ protected: 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 index 5989af74e..cfe439aed 100644 --- 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 index c0fb170f6..a15900489 100644 --- a/interfaces/cython/cantera/examples/onedim/ion_flame.py +++ b/interfaces/cython/cantera/examples/onedim/ion_flame.py @@ -19,19 +19,15 @@ gas.TPX = Tin, p, reactants # 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 index fac2c9716..ea27c8bcf 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -570,20 +570,10 @@ class IonFlame(FreeFlame): 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 @@ class IonFlame(FreeFlame): 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 index 6abbfa076..6e62d3c66 100644 --- a/interfaces/cython/cantera/onedim.pyx +++ b/interfaces/cython/cantera/onedim.pyx @@ -506,16 +506,6 @@ cdef class IonFlow(_FlowBase): else: (self.flow).fixElectricPotential() - property velocity_enabled: - """ Determines whether or not to solve the velocity.""" - def __get__(self): - return (self.flow).doVelocity(0) - def __set__(self, enable): - if enable: - (self.flow).solveVelocity() - else: - (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 index d147f2f6d..d28625601 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -946,14 +946,8 @@ class TestIonFlame(utilities.CanteraTest): # 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 index 92ec5ae06..f24f69cb7 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -48,7 +48,6 @@ 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){ @@ -56,9 +55,7 @@ void IonFlow::resize(size_t components, size_t 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(); - } } }