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

[Reactor] Independent variables for Reactor are now T,V, and Yk

Browse source 
This formulation for the reactor's governing equations significantly improves
the performance of integrator, mostly by improving the quality of the
Jacobian. The effect is small for smaller mechanisms (GRI 3.0) but can lead to
order-of-magnitude improvements for mechanisms with hundreds or thousands of
species.

Since this set of variables corresponds to the intrinsic state variables used
for IdealGasPhase objects, we also eliminate the need to iterate when setting
the state of the thermo object.

Additionally, by using temperature as an independent variable, the
temperature-dependent parts of the kinetic rate expressions do not need to be
recomputed while updating the Jacobian (this optimization is not currently
implemented).
This commit is contained in:
Ray Speth 2013-05-23 19:33:16 +00:00
parent 0ee8ac54fd
commit a59309e81e
2 changed files with 77 additions and 90 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
include/cantera/zeroD/Reactor.h
View file

@ -129,8 +129,11 @@ protected:
//! Tolerance on the temperature
doublereal m_vdot, m_Q;
doublereal m_mass; //!< total mass
vector_fp m_work;
vector_fp m_sdot; // surface production rates
vector_fp m_wdot; //!< Species net molar production rates
vector_fp m_uk; //!< Species molar internal energies
bool m_chem;
bool m_energy;
size_t m_nv;

164
src/zeroD/Reactor.cpp
View file

@ -39,30 +39,27 @@ void Reactor::getInitialConditions(double t0, size_t leny, double* y)
}
m_thermo->restoreState(m_state);
// total mass
doublereal mass = m_thermo->density() * m_vol;
// set components y + 2 ... y + K + 1 to the
// mass M_k of each species
m_thermo->getMassFractions(y+2);
scale(y + 2, y + m_nsp + 2, y + 2, mass);
// set the first component to the total internal
// energy
y[0] = m_thermo->intEnergy_mass() * mass;
// set the first component to the total mass
m_mass = m_thermo->density() * m_vol;
y[0] = m_mass;
// set the second component to the total volume
y[1] = m_vol;
// Set the third component to the temperature
y[2] = m_thermo->temperature();
// set components y+3 ... y+K+2 to the mass fractions of each species
m_thermo->getMassFractions(y+3);
// set the remaining components to the surface species
// coverages on the walls
size_t loc = m_nsp + 2;
size_t loc = m_nsp + 3;
SurfPhase* surf;
for (size_t m = 0; m < m_nwalls; m++) {
surf = m_wall[m]->surface(m_lr[m]);
if (surf) {
m_wall[m]->getCoverages(m_lr[m], y + loc);
//surf->getCoverages(y+loc);
loc += surf->nSpecies();
}
}
@ -75,7 +72,9 @@ void Reactor::initialize(doublereal t0)
{
m_thermo->restoreState(m_state);
m_sdot.resize(m_nsp, 0.0);
m_nv = m_nsp + 2;
m_wdot.resize(m_nsp, 0.0);
m_uk.resize(m_nsp, 0.0);
m_nv = m_nsp + 3;
for (size_t w = 0; w < m_nwalls; w++)
if (m_wall[w]->surface(m_lr[w])) {
m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
@ -125,44 +124,21 @@ size_t Reactor::nSensParams()
void Reactor::updateState(doublereal* y)
{
for (size_t i = 0; i < m_nsp+2; i++) {
for (size_t i = 0; i < m_nv; i++) {
AssertFinite(y[i], "Reactor::updateState",
"y[" + int2str(i) + "] is not finite");
}
// The components of y are [0] the total internal energy,
// [1] the total volume, and [2...K+2] the mass of each species.
// The components of y are [0] the total mass, [1] the total volume,
// [2] the temperature, [3...K+3] are the mass fractions of each species,
// and [K+3...] are the coverages of surface species on each wall.
m_mass = y[0];
m_vol = y[1];
// Set the mass fractions
doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0);
m_thermo->setMassFractions(y+2);
m_thermo->setMassFractions_NoNorm(y+3);
m_thermo->setState_TR(y[2], m_mass / m_vol);
if (m_energy) {
// Use Newton's method to determine the mixture temperature. Tight
// tolerances are required both for Jacobian evaluation and for
// sensitivity analysis to work correctly.
doublereal U = y[0];
doublereal T = temperature();
double dT = 100;
int i = 0;
while (abs(dT / T) > 10 * DBL_EPSILON) {
m_thermo->setState_TR(T, mass / m_vol);
double dUdT = m_thermo->cv_mass() * mass;
dT = (m_thermo->intEnergy_mass() * mass - U) / dUdT;
T -= dT;
i++;
if (i > 100) {
throw CanteraError("Reactor::updateState", "no convergence");
}
}
} else {
m_thermo->setDensity(mass/m_vol);
}
size_t loc = m_nsp + 2;
size_t loc = m_nsp + 3;
SurfPhase* surf;
for (size_t m = 0; m < m_nwalls; m++) {
surf = m_wall[m]->surface(m_lr[m]);
@ -206,9 +182,14 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
m_vdot = 0.0;
m_Q = 0.0;
double mcvdTdt = 0.0; // m * c_v * dT/dt
double dmdt = 0.0; // dm/dt (gas phase)
double* dYdt = ydot + 3;
m_thermo->getPartialMolarIntEnergies(&m_uk[0]);
// compute wall terms
size_t loc = m_nsp+2;
size_t loc = m_nsp+3;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
@ -240,68 +221,67 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
}
}
// volume equation
ydot[1] = m_vdot;
/* species equations
* Equation is:
* \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in}
* - \dot m_{out} Y_{k} + A \dot s_k.
*/
const vector_fp& mw = m_thermo->molecularWeights();
const doublereal* Y = m_thermo->massFractions();
if (m_chem) {
m_kin->getNetProductionRates(ydot+2); // "omega dot"
} else {
fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
}
for (size_t n = 0; n < m_nsp; n++) {
ydot[n+2] *= m_vol; // moles/s/m^3 -> moles/s
ydot[n+2] += m_sdot[n];
ydot[n+2] *= mw[n];
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}
/*
* Energy equation.
* \f[
* \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in}
* - \dot m_{out} h.
* \f]
*/
if (m_energy) {
ydot[0] = - m_thermo->pressure() * m_vdot - m_Q;
} else {
ydot[0] = 0.0;
double mdot_surf = 0.0; // net mass flux from surfaces
for (size_t k = 0; k < m_nsp; k++) {
// production in gas phase and from surfaces
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
mdot_surf += m_sdot[k] * mw[k];
}
dmdt += mdot_surf;
// compression work and external heat transfer
mcvdTdt += - m_pressure * m_vdot - m_Q;
for (size_t n = 0; n < m_nsp; n++) {
// heat release from gas phase and surface reations
mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol;
mcvdTdt -= m_sdot[n] * m_uk[n];
// dilution by net surface mass flux
dYdt[n] -= Y[n] * mdot_surf / m_mass;
}
// add terms for open system
if (m_open) {
const doublereal* mf = m_thermo->massFractions();
doublereal enthalpy = m_thermo->enthalpy_mass();
// outlets
for (size_t i = 0; i < m_nOutlets; i++) {
double mdot_out = m_outlet[i]->massFlowRate(time);
for (size_t n = 0; n < m_nsp; n++) {
ydot[2+n] -= mdot_out * mf[n];
}
if (m_energy) {
ydot[0] -= mdot_out * enthalpy;
}
dmdt -= mdot_out; // mass flow out of system
mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work
}
// inlets
for (size_t i = 0; i < m_nInlets; i++) {
double mdot_in = m_inlet[i]->massFlowRate(time);
dmdt += mdot_in; // mass flow into system
mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
for (size_t n = 0; n < m_nsp; n++) {
ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n);
}
if (m_energy) {
ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
// flow of species into system and dilution by other species
dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
// In combintion with h_in*mdot_in, flow work plus thermal
// energy carried with the species
mcvdTdt -= m_uk[n] / mw[n] * mdot_spec;
}
}
}
for (size_t i = 0; i < m_nsp+2; i++) {
ydot[0] = dmdt;
ydot[1] = m_vdot;
if (m_energy) {
ydot[2] = mcvdTdt / (m_mass * m_thermo->cv_mass());
} else {
ydot[2] = 0;
}
for (size_t i = 0; i < m_nv; i++) {
AssertFinite(ydot[i], "Reactor::evalEqs",
"ydot[" + int2str(i) + "] is not finite");
}
@ -350,16 +330,20 @@ std::vector<std::pair<void*, int> > Reactor::getSensitivityOrder() const
size_t Reactor::componentIndex(const string& nm) const
{
if (nm == "U") {
if (nm == "m") {
return 0;
}
if (nm == "V") {
return 1;
}
if (nm == "T") {
return 2;
}
// check for a gas species name
size_t k = m_thermo->speciesIndex(nm);
if (k != npos) {
return k + 2;
return k + 3;
}
// check for a wall species
@ -371,7 +355,7 @@ size_t Reactor::componentIndex(const string& nm) const
th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
k = th->speciesIndex(nm);
if (k != npos) {
return k + 2 + m_nsp + walloffset;
return k + 3 + m_nsp + walloffset;
} else {
walloffset += th->nSpecies();
}