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