[1D/Python] Add IonFlow to Python interface, with example and test

This commit is contained in:
bangshiuh 2017-02-09 07:51:54 -05:00 committed by Ray Speth
parent 3b12c6d662
commit e2f718c65b
8 changed files with 519 additions and 225 deletions

View file

@ -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
* methaneoxygen 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<size_t> 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<bool> 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<bool> m_do_poisson;
//! flag for solving the velocity or not
std::vector<bool> m_do_velocity;
// !electrical properties
//! electrical properties
vector_int m_speciesCharge;
// !index of species with charges
//! index of species with charges
std::vector<size_t> m_kCharge;
// !index of neutral species
//! index of neutral species
std::vector<size_t> 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

View file

@ -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

View file

@ -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)

View file

@ -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')

View file

@ -480,6 +480,43 @@ cdef class FreeFlow(_FlowBase):
self.flow = <CxxStFlow*>(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 = <CxxStFlow*>(new CxxIonFlow(gas, thermo.n_species, 2))
def set_solvingStage(self, stage):
(<CxxIonFlow*>self.flow).setSolvingStage(stage)
def set_electricPotential(self, v_inlet, v_outlet):
(<CxxIonFlow*>self.flow).setElectricPotential(v_inlet, v_outlet)
property poisson_enabled:
""" Determines whether or not to solve the energy equation."""
def __get__(self):
return (<CxxIonFlow*>self.flow).doPoisson(0)
def __set__(self, enable):
if enable:
(<CxxIonFlow*>self.flow).solvePoissonEqn()
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):
"""
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

View file

@ -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)

View file

@ -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<size_t>(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) {

85
test/data/ch4_ion.cti Normal file
View file

@ -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])