GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[1D] Change Poisson equation to first-order difference

Browse Source
This commit is contained in:
bangshiuh
2018-08-07 22:44:17 -04:00
committed by Ray Speth
parent 3ade3335de
commit ee633f3e16
8 changed files with 64 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
include/cantera/oneD/IonFlow.h
View File

@@ -34,18 +34,16 @@ 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);
virtual void resize(size_t components, size_t points);
virtual void _finalize(const double* x);
//! set to solve Poisson's equation on a point
void solvePoissonEqn(size_t j=npos);
//! set to solve electric field on a point
void solveElectricField(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];
void fixElectricField(size_t j=npos);
bool doElectricField(size_t j) {
return m_do_electric_field[j];
}
/**
@@ -66,7 +64,7 @@ public:
protected:
/*!
* This function overloads the original function. The residual function
* of Poisson's equation is added.
* of electric field is added.
*/
virtual void evalResidual(double* x, double* rsd, int* diag,
double rdt, size_t jmin, size_t jmax);
@@ -74,10 +72,11 @@ 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 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;
//! Solving phase two: the electric field equation is added coupled
//! by the electrical drift
virtual void electricFieldMethod(const double* x, size_t j0, size_t j1);
//! flag for solving electric field or not
std::vector<bool> m_do_electric_field;
//! flag for importing transport of electron
bool m_import_electron_transport;
@@ -101,33 +100,16 @@ protected:
//! solving stage
int m_stage;
//! The voltage
double m_inletVoltage;
double m_outletVoltage;
//! index of electron
size_t m_kElectron;
//! fixed electric potential value
vector_fp m_fixedElecPoten;
//! The fixed electric potential value at point j
double phi_fixed(size_t j) const {
return m_fixedElecPoten[j];
}
//! electric potential
double phi(const double* x, size_t j) const {
return x[index(c_offset_P, j)];
}
//! electric field
double E(const double* x, size_t j) const {
return -(phi(x,j+1)-phi(x,j))/(z(j+1)-z(j));
return x[index(c_offset_E, j)];
}
double dEdz(const double* x, size_t j) const {
return 2*(E(x,j)-E(x,j-1))/(z(j+1)-z(j-1));
return (E(x,j)-E(x,j-1))/(z(j)-z(j-1));
}
//! number density

2
include/cantera/oneD/StFlow.h
View File

@@ -23,7 +23,7 @@ const size_t c_offset_U = 0; // axial velocity
const size_t c_offset_V = 1; // strain rate
const size_t c_offset_T = 2; // temperature
const size_t c_offset_L = 3; // (1/r)dP/dr
const size_t c_offset_P = 4; // electric poisson's equation
const size_t c_offset_E = 4; // electric poisson's equation
const size_t c_offset_Y = 5; // mass fractions
class Transport;

7
interfaces/cython/cantera/_cantera.pxd
View File

@@ -701,10 +701,9 @@ 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 solveElectricField()
void fixElectricField()
cbool doElectricField(size_t)
cdef extern from "cantera/oneD/Sim1D.h":

33
interfaces/cython/cantera/onedim.py
View File

@@ -578,52 +578,35 @@ class IonFlameBase(FlameBase):
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)
'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] +
writer.writerow([z[n], u[n], V[n], T[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):
def electric_field_enabled(self):
""" Get/Set whether or not to solve the Poisson's equation."""
return self.flame.poisson_enabled
return self.flame.electric_field_enabled
@poisson_enabled.setter
def poisson_enabled(self, enable):
self.flame.poisson_enabled = enable
@property
def phi(self):
"""
Array containing the electric potential at each point.
"""
return self.profile(self.flame, 'ePotential')
@electric_field_enabled.setter
def electric_field_enabled(self, enable):
self.flame.electric_field_enabled = enable
@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
return self.profile(self.flame, 'eField')
def solve(self, loglevel=1, refine_grid=True, auto=False, stage=1, enable_energy=True):
self.flame.set_solvingStage(stage)

11
interfaces/cython/cantera/onedim.pyx
View File

@@ -549,18 +549,15 @@ cdef class IonFlow(_FlowBase):
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:
property electric_field_enabled:
""" Determines whether or not to solve the energy equation."""
def __get__(self):
return (<CxxIonFlow*>self.flow).doPoisson(0)
return (<CxxIonFlow*>self.flow).doElectricField(0)
def __set__(self, enable):
if enable:
(<CxxIonFlow*>self.flow).solvePoissonEqn()
(<CxxIonFlow*>self.flow).solveElectricField()
else:
(<CxxIonFlow*>self.flow).fixElectricPotential()
(<CxxIonFlow*>self.flow).fixElectricField()
cdef class Sim1D:

4
interfaces/cython/cantera/test/test_onedim.py
View File

@@ -961,7 +961,7 @@ class TestIonFreeFlame(utilities.CanteraTest):
self.sim.solve(loglevel=0, stage=2, enable_energy=True)
# Regression test
self.assertNear(max(self.sim.E), 132.1922, 1e-3)
self.assertNear(max(self.sim.E), 136.1234, 1e-3)
class TestIonBurnerFlame(utilities.CanteraTest):
@@ -986,4 +986,4 @@ class TestIonBurnerFlame(utilities.CanteraTest):
self.sim.solve(loglevel=0, stage=2, enable_energy=True)
# Regression test
self.assertNear(max(self.sim.E), 469.7287, 1e-3)
self.assertNear(max(self.sim.E), 452.9473, 1e-3)

71
src/oneD/IonFlow.cpp
View File

@@ -19,8 +19,6 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
StFlow(ph, nsp, points),
m_import_electron_transport(false),
m_stage(1),
m_inletVoltage(0.0),
m_outletVoltage(0.0),
m_kElectron(npos)
{
// make a local copy of species charge
@@ -43,23 +41,22 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
}
// no bound for electric potential
setBounds(c_offset_P, -1.0e20, 1.0e20);
setBounds(c_offset_E, -1.0e20, 1.0e20);
// Set tighter negative species limit on electron concentration to avoid
// instabilities
setBounds(c_offset_Y + m_kElectron, -1e-16, 1.0);
m_refiner->setActive(c_offset_P, false);
m_refiner->setActive(c_offset_E, false);
m_mobility.resize(m_nsp*m_points);
m_do_poisson.resize(m_points,false);
m_do_electric_field.resize(m_points,false);
}
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_fixedElecPoten.resize(m_points,0.0);
m_do_electric_field.resize(m_points,false);
}
void IonFlow::updateTransport(double* x, size_t j0, size_t j1)
@@ -83,7 +80,7 @@ void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
frozenIonMethod(x,j0,j1);
}
if (m_stage == 2) {
poissonEqnMethod(x,j0,j1);
electricFieldMethod(x,j0,j1);
}
}
@@ -114,7 +111,7 @@ void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1)
}
}
void IonFlow::poissonEqnMethod(const double* x, size_t j0, size_t j1)
void IonFlow::electricFieldMethod(const double* x, size_t j0, size_t j1)
{
for (size_t j = j0; j < j1; j++) {
double wtm = m_wtm[j];
@@ -162,17 +159,10 @@ void IonFlow::setSolvingStage(const size_t stage)
throw CanteraError("IonFlow::updateDiffFluxes",
"solution stage must be set to: "
"1) frozenIonMethod, "
"2) poissonEqnMethod");
"2) electricFieldEqnMethod");
}
}
void IonFlow::setElectricPotential(const double v1, const double v2)
{
// This method can be used when you want to add external voltage
m_inletVoltage = v1;
m_outletVoltage = v2;
}
void IonFlow::evalResidual(double* x, double* rsd, int* diag,
double rdt, size_t jmin, size_t jmax)
{
@@ -189,70 +179,70 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag,
for (size_t k : m_kCharge) {
rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1);
}
rsd[index(c_offset_P, j)] = m_inletVoltage - phi(x,j);
diag[index(c_offset_P, j)] = 0;
rsd[index(c_offset_E, j)] = E(x,0);
diag[index(c_offset_E, j)] = 0;
} else if (j == m_points - 1) {
rsd[index(c_offset_P, j)] = m_outletVoltage - phi(x,j);
diag[index(c_offset_P, j)] = 0;
rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
diag[index(c_offset_E, j)] = 0;
} else {
//-----------------------------------------------
// Poisson's equation
// Electric field by Gauss's law
//
// dE/dz = e/eps_0 * sum(q_k*n_k)
//
// E = -dV/dz
//-----------------------------------------------
rsd[index(c_offset_P, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
diag[index(c_offset_P, j)] = 0;
rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
diag[index(c_offset_E, j)] = 0;
}
}
}
void IonFlow::solvePoissonEqn(size_t j)
void IonFlow::solveElectricField(size_t j)
{
bool changed = false;
if (j == npos) {
for (size_t i = 0; i < m_points; i++) {
if (!m_do_poisson[i]) {
if (!m_do_electric_field[i]) {
changed = true;
}
m_do_poisson[i] = true;
m_do_electric_field[i] = true;
}
} else {
if (!m_do_poisson[j]) {
if (!m_do_electric_field[j]) {
changed = true;
}
m_do_poisson[j] = true;
m_do_electric_field[j] = 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);
m_refiner->setActive(c_offset_E, true);
if (changed) {
needJacUpdate();
}
}
void IonFlow::fixElectricPotential(size_t j)
void IonFlow::fixElectricField(size_t j)
{
bool changed = false;
if (j == npos) {
for (size_t i = 0; i < m_points; i++) {
if (m_do_poisson[i]) {
if (m_do_electric_field[i]) {
changed = true;
}
m_do_poisson[i] = false;
m_do_electric_field[i] = false;
}
} else {
if (m_do_poisson[j]) {
if (m_do_electric_field[j]) {
changed = true;
}
m_do_poisson[j] = false;
m_do_electric_field[j] = false;
}
m_refiner->setActive(c_offset_U, false);
m_refiner->setActive(c_offset_V, false);
m_refiner->setActive(c_offset_T, false);
m_refiner->setActive(c_offset_P, false);
m_refiner->setActive(c_offset_E, false);
if (changed) {
needJacUpdate();
}
@@ -279,14 +269,9 @@ void IonFlow::_finalize(const double* x)
{
StFlow::_finalize(x);
bool p = m_do_poisson[0];
for (size_t j = 0; j < m_points; j++) {
if (!p) {
m_fixedElecPoten[j] = phi(x, j);
}
}
bool p = m_do_electric_field[0];
if (p) {
solvePoissonEqn();
solveElectricField();
}
}

10
src/oneD/StFlow.cpp
View File

@@ -398,15 +398,15 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
// set residual of poisson's equ to zero
rsd[index(c_offset_P, 0)] = x[index(c_offset_P, j)];
rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)];
} else if (j == m_points - 1) {
evalRightBoundary(x, rsd, diag, rdt);
// set residual of poisson's equ to zero
rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)];
rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
} else { // interior points
evalContinuity(j, x, rsd, diag, rdt);
// set residual of poisson's equ to zero
rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)];
rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
//------------------------------------------------
// Radial momentum equation
@@ -582,7 +582,7 @@ string StFlow::componentName(size_t n) const
case 3:
return "lambda";
case 4:
return "ePotential";
return "eField";
default:
if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
return m_thermo->speciesName(n - c_offset_Y);
@@ -602,7 +602,7 @@ size_t StFlow::componentIndex(const std::string& name) const
return 2;
} else if (name=="lambda") {
return 3;
} else if (name == "ePotential") {
} else if (name == "eField") {
return 4;
} else {
for (size_t n=c_offset_Y; n<m_nsp+c_offset_Y; n++) {