research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[1D] Remove charge neutrality solver option from IonFlow

Browse source
This commit is contained in:
bangshiuh 2018-05-18 16:45:38 -04:00 • committed by Ray Speth
parent 3a0f46eb56
commit eeb27d84a9
7 changed files with 10 additions and 201 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
include/cantera/oneD/IonFlow.h
View file

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

3
interfaces/cython/cantera/_cantera.pxd
View file

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

6
interfaces/cython/cantera/examples/onedim/ion_flame.py
View file

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

21
interfaces/cython/cantera/onedim.py
View file

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

10
interfaces/cython/cantera/onedim.pyx
View file

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

10
interfaces/cython/cantera/test/test_onedim.py
View file

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

131
src/oneD/IonFlow.cpp
View file

@ -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();
}
}
}