From a59309e81e7218cdc0d0e747c7aa1dcdeb22bc01 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 23 May 2013 19:33:16 +0000 Subject: [PATCH] [Reactor] Independent variables for Reactor are now T,V, and Yk 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). --- include/cantera/zeroD/Reactor.h | 3 + src/zeroD/Reactor.cpp | 164 ++++++++++++++------------------ 2 files changed, 77 insertions(+), 90 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 144911070..4c077b6ff 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -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; diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index d5f32c349..91c598a01 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -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 > 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(); }