diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index 023fe38cb..75cf6b1cf 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -3,128 +3,144 @@ // This file is part of Cantera. See License.txt in the top-level directory or // at http://www.cantera.org/license.txt for license and copyright information. -#include "Domain1D.h" -#include "cantera/base/Array.h" -#include "cantera/thermo/IdealGasPhase.h" -#include "cantera/kinetics/Kinetics.h" +#ifndef CT_IONFLOW_H +#define CT_IONFLOW_H + #include "cantera/oneD/StFlow.h" -#include "cantera/oneD/Sim1D.h" -#include "cantera/IdealGasMix.h" namespace Cantera { /** - * A class for ion flow. + * This class models the ion transportation in a flame. There are three + * stages of the simulation. + * + * The first stage turns off the diffusion of ions due to the fast + * 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 + * 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. + * Reference: + * Pederson, Timothy, and R. C. Brown. + * "Simulation of electric field effects in premixed methane flames." + * Combustion and Flames 94.4(1993): 433-448. * @ingroup onedim */ class IonFlow : public FreeFlame { public: IonFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1); + //! set the solving stage + virtual void setSolvingStage(const size_t phase); + //! set electric voltage at inlet and outlet + virtual void setElectricPotential(const double v1, const double v2); - //! Turn electric field effect on/off - virtual void enableElectric(bool withElectric); - bool withElectric() const { - return m_do_electric; - } - - virtual void setSolvingPhase(const size_t phase); - - std::vector chargeList() const { - return m_kCharge; - } - - virtual void eval(size_t jg, doublereal* xg, - doublereal* rg, integer* diagg, doublereal rdt); + virtual void eval(size_t jg, double* xg, + double* rg, integer* diagg, double rdt); virtual void resize(size_t components, size_t points); - virtual void _finalize(const doublereal* x); - - void solveSpeciesEqn(size_t k=npos); - void fixSpeciesMassFrac(size_t k=npos); + virtual void _finalize(const double* x); + //! set to solve Poisson's equation on a point void solvePoissonEqn(size_t j=npos); + //! set to fix voltage on a point void fixElectricPotential(size_t j=npos); + 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]; + } protected: - virtual void updateTransport(doublereal* x, size_t j0, size_t j1); - virtual void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1); - virtual void evalPoisson(size_t j, doublereal* x, doublereal* r, integer* diag, doublereal rdt); - virtual void phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1); - virtual void phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1); - virtual void phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1); - - bool m_do_electric; - std::vector m_do_velocity; + virtual void updateTransport(double* x, size_t j0, size_t j1); + virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1); + //! evaluate the residual for Poisson's equation + virtual void evalPoisson(size_t j, double* x, double* r, integer* diag, double rdt); + //! 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; - // !electrical properties + //! electrical properties vector_int m_speciesCharge; - // !index of species with charges + //! index of species with charges std::vector m_kCharge; - // !index of neutral species + //! index of neutral species std::vector m_kNeutral; - // mobility - vector_fp m_mobi; + //! mobility + vector_fp m_mobility; - // mass fraction of ion by equlibrium - Array2D m_yCharge; + //! solving stage + int m_stage; - // IonFlow solving phase - int m_solnPhase; + //! The voltage + double m_inletVoltage; + double m_outletVoltage; - // !index of electron + //! index of electron size_t m_kElectron; - // fixed mass fraction value - vector_fp m_fixedMassFrac; - - // fixed electric potential value + //! fixed electric potential value vector_fp m_fixedElecPoten; - // fixed velocity value + //! fixed velocity value vector_fp m_fixedVelocity; //! The fixed electric potential value at point j - doublereal phi_fixed(size_t j) const { + double phi_fixed(size_t j) const { return m_fixedElecPoten[j]; } - //! The fixed mass fraction value at point j. - doublereal Y_fixed(size_t k, size_t j) const { - return m_fixedMassFrac[m_points*k+j]; + //! The fixed velocity value at point j + double u_fixed(size_t j) const { + return m_fixedVelocity[j]; } - //! The fixed velocity value at point j - doublereal u_fixed(size_t j) const { - return m_fixedVelocity[j]; - } - - // electric potential - doublereal phi(const doublereal* x, size_t j) const { + //! electric potential + double phi(const double* x, size_t j) const { return x[index(c_offset_P, j)]; } - //electric field - doublereal E(const doublereal* x, size_t j) const { + //! electric field + double E(const double* x, size_t j) const { return -(phi(x,j+1)-phi(x,j))/(z(j+1)-z(j)); } - doublereal dEdz(const doublereal* x, size_t j) const { + double dEdz(const double* x, size_t j) const { return 2*(E(x,j)-E(x,j-1))/(z(j+1)-z(j-1)); } - // number density - doublereal ND(const doublereal* x, size_t k, size_t j) const { + //! number density + double ND(const double* x, size_t k, size_t j) const { return Avogadro * m_rho[j] * Y(x,k,j) / m_wt[k]; } }; } - +#endif diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index fff940177..e6fd026e5 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -692,6 +692,19 @@ cdef extern from "cantera/oneD/StFlow.h": CxxAxiStagnFlow(CxxIdealGasPhase*, int, int) +cdef extern from "cantera/oneD/IonFlow.h": + cdef cppclass CxxIonFlow "Cantera::IonFlow": + CxxIonFlow(CxxIdealGasPhase*, int, int) + void setSolvingStage(int) + void setElectricPotential(const double, const double) + void solvePoissonEqn() + void fixElectricPotential() + cbool doPoisson(size_t) + void solveVelocity() + void fixVelocity() + cbool doVelocity(size_t) + + cdef extern from "cantera/oneD/Sim1D.h": cdef cppclass CxxSim1D "Cantera::Sim1D": CxxSim1D(vector[CxxDomain1D*]&) except +translate_exception @@ -1024,6 +1037,9 @@ cdef class _FlowBase(Domain1D): cdef class FreeFlow(_FlowBase): pass +cdef class IonFlow(_FlowBase): + pass + cdef class AxisymmetricStagnationFlow(_FlowBase): pass diff --git a/interfaces/cython/cantera/examples/onedim/ion_flame.py b/interfaces/cython/cantera/examples/onedim/ion_flame.py new file mode 100644 index 000000000..ead8d4ad4 --- /dev/null +++ b/interfaces/cython/cantera/examples/onedim/ion_flame.py @@ -0,0 +1,42 @@ +""" +A freely-propagating, premixed hydrogen flat flame with multicomponent +transport properties. +""" + +import cantera as ct +import numpy as np + +# Simulation parameters +p = ct.one_atm # pressure [Pa] +Tin = 300.0 # unburned gas temperature [K] +reactants = 'CH4:1, O2:2, N2:7.52' # premixed gas composition +width = 0.05 # m +loglevel = 1 # amount of diagnostic output (0 to 8) + +# IdealGasMix object used to compute mixture properties, set to the state of the +# upstream fuel-air mixture +gas = ct.Solution('gri30_ion.xml') +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.show_solution() + +# phase one +f.solve(loglevel=loglevel, auto=True) + +# phase two +f.solve(loglevel=loglevel, stage=2, enable_energy=False) +f.solve(loglevel=loglevel, stage=2, enable_energy=True) + +# phase 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])) + +# write the velocity, temperature, density, and mole fractions to a CSV file +f.write_csv('CH4_adiabatic.csv', quiet=False) + diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index fa57d37f6..6766af21c 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -396,7 +396,9 @@ class FreeFlame(FlameBase): """ self.inlet = Inlet1D(name='reactants', phase=gas) self.outlet = Outlet1D(name='products', phase=gas) - self.flame = FreeFlow(gas, name='flame') + if not hasattr(self, 'flame'): + # Create flame domain if not already instantiated by a child class + self.flame = FreeFlow(gas, name='flame') if width is not None: grid = np.array([0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0]) * width @@ -456,7 +458,7 @@ class FreeFlame(FlameBase): locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]]) def get_flame_speed_reaction_sensitivities(self): - r""" + """ Compute the normalized sensitivities of the laminar flame speed :math:`S_u` with respect to the reaction rate constants :math:`k_i`: @@ -484,6 +486,103 @@ class FreeFlame(FlameBase): return self.solve_adjoint(perturb, self.gas.n_reactions, dgdx) / Su0 +class IonFlame(FreeFlame): + __slots__ = ('inlet', 'outlet', 'flame') + + def __init__(self, gas, grid=None, width=None): + self.flame = IonFlow(gas, name='flame') + 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 + 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) + + def write_csv(self, filename, species='X', quiet=True): + """ + Write the velocity, temperature, density, electric potential, + , electric field stregth, and species profiles to a CSV file. + + :param filename: + Output file name + :param species: + Attribute to use obtaining species profiles, e.g. ``X`` for + mole fractions or ``Y`` for mass fractions. + """ + z = self.grid + T = self.T + u = self.u + V = self.V + phi = self.phi + E = self.E + + csvfile = open(filename, 'w') + writer = _csv.writer(csvfile) + writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', 'T (K)', + 'phi (V)', 'E (V/m)', 'rho (kg/m3)'] + self.gas.species_names) + for n in range(self.flame.n_points): + self.set_gas_state(n) + writer.writerow([z[n], u[n], V[n], T[n], phi[n], E[n], self.gas.density] + + list(getattr(self.gas, species))) + csvfile.close() + if not quiet: + print("Solution saved to '{0}'.".format(filename)) + + @property + def poisson_enabled(self): + """ Get/Set whether or not to solve the energy equation.""" + return self.flame.poisson_enabled + + @poisson_enabled.setter + def poisson_enabled(self, enable): + self.flame.poisson_enabled = enable + + @property + def velocity_enabled(self): + """ Get/Set whether or not to solve the energy equation.""" + return self.flame.velocity_enabled + + @velocity_enabled.setter + def velocity_enabled(self, enable): + self.flame.velocity_enabled = enable + + @property + def phi(self): + """ + Array containing the electric potential at each point. + """ + return self.profile(self.flame, 'ePotential') + + @property + def E(self): + """ + Array containing the electric field strength at each point. + """ + z = self.grid + phi = self.phi + np = self.flame.n_points + Efield = [] + Efield.append((phi[0] - phi[1]) / (z[1] - z[0])) + # calculate E field strength + for n in range(1,np-1): + Efield.append((phi[n-1] - phi[n+1]) / (z[n+1] - z[n-1])) + Efield.append((phi[np-2] - phi[np-1]) / (z[np-1] - z[np-2])) + return Efield + + class BurnerFlame(FlameBase): """A burner-stabilized flat flame.""" __slots__ = ('burner', 'flame', 'outlet') diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx index 3388ed47b..b0dc26a43 100644 --- a/interfaces/cython/cantera/onedim.pyx +++ b/interfaces/cython/cantera/onedim.pyx @@ -480,6 +480,43 @@ cdef class FreeFlow(_FlowBase): self.flow = (new CxxFreeFlame(gas, thermo.n_species, 2)) +cdef class IonFlow(_FlowBase): + """ + An ion flow domain. + + In an ion flow dommain, the electric drift is added to the diffusion flux + """ + def __cinit__(self, _SolutionBase thermo, *args, **kwargs): + gas = getIdealGasPhase(thermo) + self.flow = (new CxxIonFlow(gas, thermo.n_species, 2)) + + def set_solvingStage(self, stage): + (self.flow).setSolvingStage(stage) + + def set_electricPotential(self, v_inlet, v_outlet): + (self.flow).setElectricPotential(v_inlet, v_outlet) + + property poisson_enabled: + """ Determines whether or not to solve the energy equation.""" + def __get__(self): + return (self.flow).doPoisson(0) + def __set__(self, enable): + if enable: + (self.flow).solvePoissonEqn() + 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): """ An axisymmetric flow domain. @@ -1072,7 +1109,7 @@ cdef class Sim1D: self.sim.clearStats() def solve_adjoint(self, perturb, n_params, dgdx, g=None, dp=1e-5): - r""" + """ Find the sensitivities of an objective function using an adjoint method. For an objective function :math:`g(x, p)` where :math:`x` is the state diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index 945d8f62f..20caba67b 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -853,3 +853,40 @@ class TestTwinFlame(utilities.CanteraTest): def test_case1(self): self.solve(phi=0.4, T=300, width=0.05, P=0.1) + + +class TestIonFlame(utilities.CanteraTest): + def test_ion_profile(self): + reactants = 'CH4:0.216, O2:2' + p = ct.one_atm + Tin = 300 + width = 0.03 + + # IdealGasMix object used to compute mixture properties + self.gas = ct.Solution('ch4_ion.cti') + self.gas.TPX = Tin, p, reactants + self.sim = ct.IonFlame(self.gas, width=width) + self.sim.set_refine_criteria(ratio=4, slope=0.8, curve=1.0) + + # stage one + self.sim.solve(loglevel=0, auto=True) + T1 = self.sim.T[-1] + + # stage two + self.sim.solve(loglevel=0, stage=2, enable_energy=False) + + # stage two + self.sim.solve(loglevel=0, stage=2, enable_energy=True) + Electron2 = self.sim.value(self.sim.flame, 'E', self.sim.flame.n_points-1) + + #stage three + self.sim.solve(loglevel=0, stage=3, enable_energy=True) + Electron3 = self.sim.value(self.sim.flame, 'E', self.sim.flame.n_points-1) + T3 = self.sim.T[-1] + + # check Temperature at outlet + self.assertNear(T1, T3, 1e-3) + self.assertNotEqual(T1, T3) + # check Electron concentration at outlet + self.assertNear(Electron2, Electron3, 1e-13) + self.assertNotEqual(Electron2, Electron3) diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index 226d71c97..0f7e50388 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -8,8 +8,6 @@ #include "cantera/base/ctml.h" #include "cantera/transport/TransportBase.h" #include "cantera/numerics/funcs.h" -#include "cantera/oneD/Domain1D.h" - using namespace std; @@ -18,8 +16,10 @@ namespace Cantera IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : FreeFlame(ph, nsp, points), - m_do_electric(false), - m_solnPhase(1) + m_stage(1), + m_inletVoltage(0.0), + m_outletVoltage(0.0), + m_kElectron(npos) { // make a local copy of species charge for (size_t k = 0; k < m_nsp; k++) { @@ -35,56 +35,72 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : } } - // Find the index of electron - if (m_thermo->speciesIndex("E") < m_nsp ) { + // Find the index of electron + if (m_thermo->speciesIndex("E") != npos ) { m_kElectron = m_thermo->speciesIndex("E"); + setTransientTolerances(1.0e-5, 1.0e-18, c_offset_Y + m_kElectron); + setSteadyTolerances(1.0e-5, 1.0e-16, c_offset_Y + m_kElectron); + } + if (m_thermo->speciesIndex("HCO+") != npos ) { + size_t k = m_thermo->speciesIndex("HCO+"); + setTransientTolerances(1.0e-5, 1.0e-18, c_offset_Y + k); + setSteadyTolerances(1.0e-5, 1.0e-16, c_offset_Y + k); + } + if (m_thermo->speciesIndex("H3O+") != npos ) { + size_t k = m_thermo->speciesIndex("H3O+"); + setTransientTolerances(1.0e-5, 1.0e-15, c_offset_Y + k); + setSteadyTolerances(1.0e-5, 1.0e-13, c_offset_Y + k); } // mass fraction bounds (strict bound for ions) for (size_t k : m_kCharge) { - setBounds(c_offset_Y+k, -1.0e-20, 1e-5); + setBounds(c_offset_Y+k, -1.0e-20, 1.0e5); } + // no bound for electric potential setBounds(c_offset_P, -1.0e20, 1.0e20); + m_refiner->setActive(c_offset_P, false); - - m_mobi.resize(m_nsp*m_points); + 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_mobi.resize(m_nsp*m_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_fixedMassFrac.resize(m_points*m_nsp); m_fixedElecPoten.resize(m_points,0.0); m_fixedVelocity.resize(m_points); } -void IonFlow::updateTransport(doublereal* x, size_t j0, size_t j1) +void IonFlow::updateTransport(double* x, size_t j0, size_t j1) { StFlow::updateTransport(x,j0,j1); for (size_t j = j0; j < j1; j++) { setGasAtMidpoint(x,j); - m_trans->getMobilities(&m_mobi[j*m_nsp]); - m_mobi[m_kElectron+m_nsp*j] = 0.4; - m_diff[m_kElectron+m_nsp*j] = 0.4*(Boltzmann * T(x,j)) / ElectronCharge; + m_trans->getMobilities(&m_mobility[j*m_nsp]); + if (m_kElectron != npos) { + m_mobility[m_kElectron+m_nsp*j] = 0.4; + m_diff[m_kElectron+m_nsp*j] = 0.4*(Boltzmann * T(x,j)) / ElectronCharge; + } } } -void IonFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1) +void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1) { - if (m_solnPhase == 1) { - phaseOneDiffFluxes(x,j0,j1); - } else if (m_solnPhase == 2) { - phaseTwoDiffFluxes(x,j0,j1); - } else { - phaseThreeDiffFluxes(x,j0,j1); + if (m_stage == 1) { + frozenIonMethod(x,j0,j1); + } + if (m_stage == 2) { + chargeNeutralityModel(x,j0,j1); + } + if (m_stage == 3) { + poissonEqnMethod(x,j0,j1); } } -void IonFlow::phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1) +void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1) { for (size_t j = j0; j < j1; j++) { double wtm = m_wtm[j]; @@ -98,29 +114,30 @@ void IonFlow::phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1) } // correction flux to insure that \sum_k Y_k V_k = 0. - for (size_t k : m_kNeutral) { + for (size_t k : m_kNeutral) { m_flux(k,j) += sum*Y(x,k,j); } // flux for ions + // Set flux to zero to prevent some fast charged species (e.g. electron) + // to run away for (size_t k : m_kCharge) { m_flux(k,j) = 0; } } } -void IonFlow::phaseTwoDiffFluxes(const doublereal* 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 - double sum_flux = 0.0; - for (size_t k = 0; k < m_nsp; k++) { + 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; - sum_flux -= m_flux(k,j); } // ambipolar diffusion @@ -130,83 +147,106 @@ void IonFlow::phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1) 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); - sum += m_mobi[k+m_nsp*j] * Xav * q_k * q_k; + // 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; } - double drift; - double sum_drift = 0.0; 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_mobi[k+m_nsp*j] * Xav / sum; + 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; - sum_drift -= drift; } // correction flux - for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,j) += Y(x,k,j) * sum_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::phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1) +void IonFlow::poissonEqnMethod(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 double sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { + 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; sum -= m_flux(k,j); } - // correction flux - for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,j) += Y(x,k,j) * sum; - } - // ambipolar diffusion - double drift; double E_ambi = E(x,j); - sum = 0.0; for (size_t k : m_kCharge) { double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1)); - drift = rho * Yav * E_ambi; - drift *= m_speciesCharge[k] * m_mobi[k+m_nsp*j]; + double drift = rho * Yav * E_ambi + * m_speciesCharge[k] * m_mobility[k+m_nsp*j]; m_flux(k,j) += drift; - sum -= drift; } - // correction 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) { - m_flux(k,j) += Y(x,k,j) * sum; + 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::enableElectric(bool withElectric) +void IonFlow::setSolvingStage(const size_t stage) { - m_do_electric = withElectric; + if (stage == 1 || stage == 2 || stage == 3) { + m_stage = stage; + } else { + throw CanteraError("IonFlow::updateDiffFluxes", + "solution phase must be set to:" + "1: frozenIonMethod" + "2: chargeNeutralityModel" + "3: poissonEqnMethod"); + } } -void IonFlow::setSolvingPhase(const size_t phase) +void IonFlow::setElectricPotential(const double v1, const double v2) { - m_solnPhase = phase; + // This method can be used when you want to add external voltage + m_inletVoltage = v1; + m_outletVoltage = v2; } -void IonFlow::eval(size_t jg, doublereal* xg, - doublereal* rg, integer* diagg, doublereal rdt) +void IonFlow::eval(size_t jg, double* xg, + double* rg, integer* diagg, double rdt) { StFlow::eval(jg, xg, rg, diagg, rdt); + if (m_stage != 3) { + return; + } // start of local part of global arrays - doublereal* x = xg + loc(); - doublereal* rsd = rg + loc(); + double* x = xg + loc(); + double* rsd = rg + loc(); integer* diag = diagg + loc(); - size_t jmin, jmax; if (jg == npos) { // evaluate all points jmin = 0; @@ -216,66 +256,40 @@ void IonFlow::eval(size_t jg, doublereal* xg, jmin = std::max(jpt, 1) - 1; jmax = std::min(jpt+1,m_points-1); } - // the boundary points are not applied + for (size_t j = jmin; j <= jmax; j++) { if (j == 0) { - rsd[index(c_offset_P, j)] = -phi(x,j); + rsd[index(c_offset_P, j)] = m_inletVoltage - phi(x,j); diag[index(c_offset_P, j)] = 0; - for ( size_t k : m_kCharge) { - rsd[index(c_offset_Y + k, j)] = Y(x,k,j); - diag[index(c_offset_Y + k, j)] = 0; + // set ions boundary for better convergence + for (size_t k : m_kCharge) { + rsd[index(c_offset_Y + k, j)] = Y(x,k,j+1) - Y(x,k,j); } } else if (j == m_points - 1) { - rsd[index(c_offset_P, j)] = -phi(x,j); + rsd[index(c_offset_P, j)] = m_outletVoltage - phi(x,j); diag[index(c_offset_P, j)] = 0; - for ( size_t k : m_kCharge) { - rsd[index(c_offset_Y + k, j)] = Y(x,k,j); - diag[index(c_offset_Y + k, j)] = 0; - } } else { + evalPoisson(j,x,rsd,diag,rdt); if (!m_do_velocity[j]) { + // This method is used when you disable energy equation + // but still maintain the velocity profile rsd[index(c_offset_U, j)] = u(x,j) - u_fixed(j); diag[index(c_offset_U, j)] = 0; } - for (size_t k = 0; k < m_nsp; k++) { - if (!m_do_species[k]) { - rsd[index(c_offset_Y + k, j)] = Y(x,k,j) - Y_fixed(k,j); - rsd[index(c_offset_Y + k, j)] -= rdt*(Y(x,k,j) - Y_prev(k,j)); - diag[index(c_offset_Y + k, j)] = 1; - } - } - } - } - - // convinent method due to interference - for (size_t j = jmin; j <= jmax; j++) { - if (j == 0) { - rsd[index(c_offset_P, j)] = -phi(x,j); - diag[index(c_offset_P, j)] = 0; - } else if (j == m_points - 1) { - rsd[index(c_offset_P, j)] = -phi(x,j); - diag[index(c_offset_P, j)] = 0; - } else { - if (m_do_poisson[j]) { - evalPoisson(j,x,rsd,diag,rdt); - } else { - rsd[index(c_offset_P, j)] = phi(x,j) - phi_fixed(j); - diag[index(c_offset_P, j)] = 0; - } } } } -void IonFlow::evalPoisson(size_t j, doublereal* x, doublereal* rsd, integer* diag, doublereal rdt) +void IonFlow::evalPoisson(size_t j, double* x, double* rsd, integer* diag, double rdt) { //----------------------------------------------- // Poisson's equation // // dE/dz = e/eps_0 * sum(q_k*n_k) // - // E = -dV/dz + // E = -dV/dz //----------------------------------------------- - doublereal chargeDensity = 0.0; + double chargeDensity = 0.0; for (size_t k : m_kCharge) { chargeDensity += m_speciesCharge[k] * ElectronCharge * ND(x,k,j); } @@ -283,48 +297,6 @@ void IonFlow::evalPoisson(size_t j, doublereal* x, doublereal* rsd, integer* dia diag[index(c_offset_P, j)] = 0; } -void IonFlow::solveSpeciesEqn(size_t k) -{ - bool changed = false; - if (k == npos) { - for (size_t i = 0; i < m_nsp; i++) { - if (!m_do_energy[i]) { - changed = true; - } - m_do_species[i] = true; - } - } else { - if (!m_do_species[k]) { - changed = true; - } - m_do_species[k] = true; - } - if (changed) { - needJacUpdate(); - } -} - -void IonFlow::fixSpeciesMassFrac(size_t k) -{ - bool changed = false; - if (k == npos) { - for (size_t i = 0; i < m_nsp; i++) { - if (m_do_species[i]) { - changed = true; - } - m_do_species[i] = false; - } - } else { - if (m_do_species[k]) { - changed = true; - } - m_do_species[k] = false; - } - if (changed) { - needJacUpdate(); - } -} - void IonFlow::solvePoissonEqn(size_t j) { bool changed = false; @@ -341,9 +313,10 @@ void IonFlow::solvePoissonEqn(size_t j) } m_do_poisson[j] = true; } - m_refiner->setActive(0, true); - m_refiner->setActive(1, true); - m_refiner->setActive(2, true); + m_refiner->setActive(c_offset_U, true); + m_refiner->setActive(c_offset_V, true); + m_refiner->setActive(c_offset_T, true); + m_refiner->setActive(c_offset_P, true); if (changed) { needJacUpdate(); } @@ -368,6 +341,7 @@ void IonFlow::fixElectricPotential(size_t j) m_refiner->setActive(0, false); m_refiner->setActive(1, false); m_refiner->setActive(2, false); + m_refiner->setActive(4, false); if (changed) { needJacUpdate(); } @@ -421,22 +395,10 @@ void IonFlow::fixVelocity(size_t j) } } -void IonFlow::_finalize(const doublereal* x) +void IonFlow::_finalize(const double* x) { FreeFlame::_finalize(x); - for (size_t k = 0; k < m_nsp; k++) { - bool y = m_do_species[k]; - if (!y) { - for (size_t j = 0; j < m_points; j++) { - m_fixedMassFrac[m_points*k+j] = Y(x,k,j); - } - } - } - - // This method is still not tested - // not sure why you want to return to original state - // if not doing on point zero bool p = m_do_poisson[0]; for (size_t j = 0; j < m_points; j++) { if (!p) { @@ -446,7 +408,7 @@ void IonFlow::_finalize(const doublereal* 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) { diff --git a/test/data/ch4_ion.cti b/test/data/ch4_ion.cti new file mode 100644 index 000000000..8b42063df --- /dev/null +++ b/test/data/ch4_ion.cti @@ -0,0 +1,85 @@ +units(length='cm', time='s', quantity='mol', act_energy='cal/mol') + +ideal_gas(name='gas', + elements='O H C N E', + species=['''gri30: H2 H O O2 OH H2O HO2 + H2O2 C CH CH2 CH2(S) CH3 CH4 + CO CO2 HCO CH2O CH3O N2''', + 'HCO+ H3O+ E'], + reactions=['gri30: all', 'all'], + transport='Mix', + options=['skip_undeclared_species', 'skip_undeclared_third_bodies'], + initial_state=state(temperature=300.0, pressure=OneAtm)) + +#------------------------------------------------------------------------------- +# Species data +#------------------------------------------------------------------------------- + +species(name = 'HCO+', + atoms = ' H:1 C:1 O:1 E:-1 ', + thermo = ( + NASA( [ 300.00, 1000.00], [ 2.473973600E+00, 8.671559000E-03, + -1.003150000E-05, 6.717052700E-09, -1.787267400E-12, + 9.914660800E+04, 8.175711870E+00] ), + NASA( [ 1000.00, 5000.00], [ 3.741188000E+00, 3.344151700E-03, + -1.239712100E-06, 2.118938800E-10, -1.370415000E-14, + 9.888407800E+04, 2.078613570E+00] ) + ), + transport=gas_transport(geom='linear', + diam=3.59, + well_depth=498.0, + polar=2.5, + rot_relax=0.0), + note = 'J12/70') + +species(name = 'H3O+', + atoms = ' H:3 O:1 E:-1 ', + thermo = ( + NASA( [ 298.15, 1000.00], [ 3.792952700E+00, -9.108540000E-04, + 1.163635490E-05, -1.213648870E-08, 4.261596630E-12, + 7.075124010E+04, 1.471568560E+00] ), + NASA( [ 1000.00, 6000.00], [ 2.496477160E+00, 5.728449200E-03, + -1.839532810E-06, 2.735774390E-10, -1.540939850E-14, + 7.097291130E+04, 7.458507790E+00] ) + ), + transport=gas_transport(geom='nonlinear', + diam=2.605, + well_depth=572.4, + dipole=1.844, + polar=1.5, + rot_relax=2.1), + note = 'TPIS89') + +species(name = 'E', + atoms = ' E:1 ', + thermo = ( + NASA( [ 200.00, 1000.00], [ 2.500000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + -7.453750000E+02, -1.172469020E+01] ), + NASA( [ 1000.00, 6000.00], [ 2.500000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + -7.453750000E+02, -1.172469020E+01] ) + ), + transport=gas_transport(geom='atom', + diam=2.05, + well_depth=145.0, + polar=0.667, + rot_relax=0.0), + note = 'gas L10/92') + +#------------------------------------------------------------------------------- +# Reaction data +#------------------------------------------------------------------------------- + +reaction('CH + O => HCO+ + E', [2.51E+11, 0.0, 1700]) + +reaction('HCO+ + H2O => H3O+ + CO', [1.51E+15, 0.0, 0.0]) + +reaction('H3O+ + E => H2O + H', [2.29E+18, -0.5, 0.0]) + +reaction('H3O+ + E => OH + H + H', [7.95E+21, -1.4, 0.0]) + +reaction('H3O+ + E => H2 + OH', [1.25E+19, -0.5, 0.0]) + +reaction('H3O+ + E => O + H2 + H', [6.0E+17, -0.3, 0.0]) +