ReactorRHS.cpp

Go to the documentation of this file.

#include "ReactorRHS.h"
#include "units.h"
 
#include "cantera/thermo/PlasmaPhase.h"
#include "cantera/kinetics/Kinetics.h"
#include "cantera/kinetics/Reaction.h"
#include "cantera/base/ct_defs.h"
 
#include <algorithm>
#include <cmath>
#include <sstream>
#include <string>
 
namespace rizer {
 
namespace {
// Exact twin of rizer.kinetics.electron_reactions.reaction_involves_electron: the
// electron participates iff it appears in the reaction *equation*. Cantera cancels a
// spectator electron (one on both sides, e.g. "CH4 + e- => CH3 + H + e-") from the
// reactant/product stoichiometry, so a species-map test would miss such electron-
// impact excitation/dissociation reactions. Both implementations replace the
// reaction-arrow and stoichiometric-plus characters '+','=','<','>' with spaces and
// then look for the electron as a standalone token; they MUST stay byte-for-byte
// equivalent so the native and Python reactors select the same electron-impact
// reactions for the inelastic electron-to-gas power. Change them together.
bool equationHasElectron(const std::string& equation, const std::string& electron)
{
    std::string s = equation;
    for (char& ch : s) {
        if (ch == '+' || ch == '=' || ch == '<' || ch == '>') ch = ' ';
    }
    std::istringstream iss(s);
    std::string token;
    while (iss >> token) {
        if (token == electron) return true;
    }
    return false;
}
} // namespace
 
// Physical constants from units.h (mirror rizer/misc/units.py) rather than
// Cantera's CODATA-2018 values, to match the Python reactor reference.
constexpr double GasConstant = units::R_kmol;  // J/kmol/K
constexpr double Boltzmann   = units::k_b;     // J/K
 
ReactorRHS::ReactorRHS(std::shared_ptr<Cantera::Solution> sol, Config cfg)
    : m_sol(std::move(sol))
    , m_nsp(m_sol->thermo()->nSpecies())
    , m_ie(m_sol->thermo()->speciesIndex("e-"))
    , m_reacting(cfg.reacting)
    , m_collision(std::move(cfg.collision))
    // The legacy fallback members (m_sigma, m_nu_m_tab, m_nu_m, m_nu_E) are always
    // built from cfg here, but transport() only reads them when m_collision.empty()
    // is true -- i.e. only one of the two closures is ever actually exercised at
    // evaluation time.
    , m_sigma(PropertyTable::fromTableOrScalar(cfg.sigma_T, cfg.sigma_v))
    , m_nu_m_tab(PropertyTable::fromTableOrScalar(cfg.nu_m_T, cfg.nu_m_v))
    , m_nu_m(cfg.nu_m)
    , m_nu_E(cfg.nu_E)
{
    if (m_ie == Cantera::npos) {
        throw std::invalid_argument("ReactorRHS: mechanism has no 'e-' species.");
    }
    m_wt.assign(m_sol->thermo()->molecularWeights().begin(),
                m_sol->thermo()->molecularWeights().end());
    m_Yc.resize(m_nsp);
    m_wdot.resize(m_nsp);
    m_uk.resize(m_nsp);
    m_cpR.resize(m_nsp);
    m_cp_mass.resize(m_nsp);
    m_X.resize(m_nsp);
    m_nk.resize(m_nsp);
    const std::size_t nR = m_sol->kinetics()->nReactions();
    m_rop.resize(nR);
    m_is_electron_rxn.assign(nR, 0);
    // Select electron-impact reactions exactly as the Python reference does: the
    // electron participates iff it appears in the reaction equation. This covers
    // two-temperature-plasma, three-body-two-temperature-plasma,
    // reverse-two-temperature-plasma, Druyvesteyn and janev-* reactions -- including
    // those with a spectator electron that Cantera cancels from the stoichiometry.
    const std::string e_name = m_sol->thermo()->speciesName(m_ie);
    bool any_electron_rxn = false;
    for (std::size_t i = 0; i < nR; i++) {
        if (equationHasElectron(m_sol->kinetics()->reaction(i)->equation(), e_name)) {
            m_is_electron_rxn[i] = 1;
            any_electron_rxn = true;
        }
    }
 
    // Fixed per-reaction threshold energy, epsilon_th_i = sum_k(nu_ki * H_k(0K)),
    // computed once here from the caller-supplied per-species 0K formation
    // enthalpies (mirrors rizer.kinetics.electron_reactions.ElectronicReactionEnergetics
    // exactly -- see Config::species_h0k). Required whenever there is at least one
    // electron-impact reaction and chemistry is active, since P_inel would otherwise
    // silently be zero.
    m_eps_th.assign(nR, 0.0);
    if (m_reacting && any_electron_rxn) {
        if (cfg.species_h0k.empty()) {
            throw std::invalid_argument(
                "ReactorRHS: Config.species_h0k is required (mechanism has an "
                "electron-impact reaction and reacting=true).");
        }
        if (cfg.species_h0k.size() != m_nsp) {
            throw std::invalid_argument(
                "ReactorRHS: Config.species_h0k must have length nSpecies().");
        }
        m_sol->kinetics()->getReactionDelta(
            Cantera::span<const double>(cfg.species_h0k.data(), m_nsp),
            Cantera::span<double>(m_eps_th.data(), m_eps_th.size()));
    }
}
 
std::string ReactorRHS::speciesName(std::size_t k) const
{
    return m_sol->thermo()->speciesName(k);
}
 
void ReactorRHS::setState(double Tg, double Te, const double* Y, double rho, double& P,
                          double& Xe, double& Tmean, double& n_e) const
{
    auto thermo = m_sol->thermo();
    // Clip any negative mass fractions (roundoff from the ODE integrator) to zero
    // before handing them to Cantera; setMassFractions_NoNorm keeps them
    // un-renormalized so Y still sums to (numerically) 1 as the integrator expects.
    for (std::size_t k = 0; k < m_nsp; k++) m_Yc[k] = (Y[k] > 0.0) ? Y[k] : 0.0;
    const double Tg_c = std::max(Tg, 200.0);
    const double Te_c = std::max(Te, 200.0);
    thermo->setMassFractions_NoNorm(Cantera::span<const double>(m_Yc.data(), m_nsp));
    thermo->setElectronTemperature(Te_c);
    thermo->setState_TD(Tg_c, rho);
    P = thermo->pressure();
    Xe = thermo->moleFraction(m_ie);
    // Two-temperature EOS mean temperature (PlasmaPhase convention),
    // @f$P = n_\text{tot} k_B T_\text{mean}@f$, with Tmean the mole-fraction-weighted
    // average of Tg and Te. n_e below uses Tmean (not Tg) to stay consistent with how
    // Cantera's PlasmaPhase computes P internally.
    Tmean = (1.0 - Xe) * Tg_c + Xe * Te_c;
    n_e = Xe * P / (Boltzmann * Tmean);
}
 
void ReactorRHS::transport(double Te, double n_e, double P, double Tmean, double& sigma,
                           double& nu_E, double& nu_eH, double& nu_eI,
                           double& nu_eH_mass_weighted) const
{
    const double Te_c = std::max(Te, 200.0);
    // Closure selection: the composition-resolved CollisionModel is preferred and used
    // exclusively whenever it was supplied (m_collision non-empty); the legacy
    // tabulated/scalar closure below is only reached when no CollisionModel was given.
    // Exactly one branch executes per call -- see the Config comment in ReactorRHS.h.
    if (!m_collision.empty()) {
        // Composition-resolved: convert Cantera's mole fractions X_k [-] to
        // per-species number densities n_k [1/m^3] via the ideal-gas number density
        // @f$n_\text{tot} = P/(k_B T_\text{mean})@f$, since CollisionModel::evaluate()
        // needs absolute densities, not fractions.
        const double n_tot = P / (Boltzmann * Tmean);
        m_sol->thermo()->getMoleFractions(Cantera::span<double>(m_X.data(), m_nsp));
        for (std::size_t k = 0; k < m_nsp; k++) m_nk[k] = m_X[k] * n_tot;
        m_collision.evaluate(Te_c, n_e, m_nk.data(), m_nsp, sigma, nu_E, nu_eH, nu_eI,
                              nu_eH_mass_weighted);
    } else {
        // Legacy closure: prefer nu_m(Te) (table or scalar) and derive sigma from the
        // Drude relation @f$\sigma = n_e e^2 / (m_e \nu_m)@f$; if no nu_m is available
        // at all (nu_m_eff <= 0), fall back further to a directly tabulated sigma(Te).
        // nu_E is always the fixed scalar in this legacy path (no Te-dependent table
        // for nu_E). nu_eH/nu_eI/nu_eH_mass_weighted need per-species composition,
        // unavailable here.
        const double nu_m_eff = m_nu_m_tab.empty() ? m_nu_m : m_nu_m_tab.eval(Te_c);
        sigma = (nu_m_eff > 0.0)
            ? n_e * units::e * units::e / (units::m_e * nu_m_eff)
            : m_sigma.eval(Te_c);
        nu_E = m_nu_E;
        nu_eH = 0.0;
        nu_eI = 0.0;
        nu_eH_mass_weighted = 0.0;
    }
}
 
void ReactorRHS::rates(double Tg, double Te, const double* Y, double rho, double E,
                       double& dTg, double& dTe, double* dY, double& lhs_Tg,
                       double& lhs_Te, Diagnostics* diag) const
{
    auto thermo = m_sol->thermo();
    auto kin = m_sol->kinetics();
 
    double P, Xe, Tmean, n_e;
    setState(Tg, Te, Y, rho, P, Xe, Tmean, n_e);
    const double Tg_c = std::max(Tg, 200.0);
    const double Te_c = std::max(Te, 200.0);
 
    thermo->getCp_R(m_cpR);
    const double We = m_wt[m_ie];
    const double cve = 1.5 * GasConstant / We;
 
    // Mass-basis species heat capacity c_p,k = (Cp_k/R)*R/W_k [J/kg/K],
    // exposed via speciesCp() for PlasmaChannel1D's enthalpy-diffusion
    // energy-equation term. The electron slot instead mirrors the analytic
    // monatomic-ideal-gas value c_p,e = cv_e + R/W_e (derived from cve just
    // above, not from getCp_R()), for consistency with how cve/lhs_Te
    // already bypass Cantera's own electron Cp elsewhere in this function.
    for (std::size_t k = 0; k < m_nsp; k++) {
        m_cp_mass[k] = m_cpR[k] * GasConstant / m_wt[k];
    }
    m_cp_mass[m_ie] = cve + GasConstant / We;
 
    // Gas-side heat capacity at constant volume, @f$c_v/R = c_p/R - 1@f$ per species
    // (ideal gas), mass-weighted over the heavy species only (electron excluded -- it
    // has its own energy equation with capacitance lhs_Te below) and converted to a
    // per-volume capacitance via rho: cv_heavy has units J/(kg K) before the rho
    // multiply.
    double cv_heavy = 0.0;
    for (std::size_t k = 0; k < m_nsp; k++) {
        if (k == m_ie) continue;
        cv_heavy += m_Yc[k] * (m_cpR[k] - 1.0) * GasConstant / m_wt[k];
    }
    // Electron-energy capacitance @f$Y_e \rho c_{v,e}@f$; Y_e is floored away from
    // zero so lhs_Te stays finite/positive even when the electron mass fraction
    // underflows to ~0.
    const double Ye_cap = std::max(m_Yc[m_ie], 1.0e-12);
    lhs_Te = Ye_cap * cve * rho;
    lhs_Tg = std::max(cv_heavy * rho, 1.0e-3);
 
    double sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted;
    transport(Te_c, n_e, P, Tmean, sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted);
 
    // Elastic electron-heavy energy exchange (classic elastic-collision power
    // transfer rate from electrons to the heavy gas, positive when Te > Tg):
    // @f[
    //     P_\text{el} = \nu_E n_e \frac{3}{2} k_B (T_e - T_g)
    // @f]
    // It appears with a minus sign in the electron balance and a plus sign in the
    // gas balance below (energy conserved between the two) [W/m^3].
    const double P_el = nu_E * n_e * 1.5 * Boltzmann * (Te_c - Tg_c);
    // Joule heating, @f$P_\text{joule} = \sigma E^2@f$ [W/m^3]: the resistive power
    // deposited into the electron population by the electric field; it only enters
    // the electron-energy balance (dTe), never the gas balance (dTg) directly.
    const double P_joule = sigma * E * E;
 
    double P_chem = 0.0, P_chem_e = 0.0, P_inel = 0.0;
    for (std::size_t k = 0; k < m_nsp; k++) dY[k] = 0.0;
    if (m_reacting) {
        kin->getNetProductionRates(m_wdot);
        kin->getNetRatesOfProgress(m_rop);
        thermo->getPartialMolarIntEnergies(m_uk);
        for (std::size_t k = 0; k < m_nsp; k++) {
            // m_wdot[k] is the net production rate in [kmol/m^3/s] (Cantera
            // convention); multiplying by the molecular weight [kg/kmol] and
            // dividing by rho [kg/m^3] converts it to the mass-fraction time
            // derivative @f$dY_k = \dot\omega_k W_k / \rho@f$ [1/s].
            dY[k] = m_wdot[k] * m_wt[k] / rho;
            if (k != m_ie) P_chem += m_uk[k] * m_wdot[k];
        }
        P_chem_e = m_wdot[m_ie] * cve * We * Te_c;
        // Fixed threshold energy per reaction (m_eps_th, precomputed once in the
        // constructor from 0K formation enthalpies -- see plasma_power_inelastic
        // in the Python reference).
        for (std::size_t i = 0; i < m_rop.size(); i++) {
            if (m_is_electron_rxn[i]) P_inel += m_eps_th[i] * m_rop[i];
        }
    }
 
    dTe = (-P_chem_e - P_el - P_inel + P_joule) / lhs_Te;
    dTg = (-P_chem + P_el + P_inel) / lhs_Tg;
 
    if (!std::isfinite(dTe)) dTe = 0.0;
    if (!std::isfinite(dTg)) dTg = 0.0;
    for (std::size_t k = 0; k < m_nsp; k++) {
        if (!std::isfinite(dY[k])) dY[k] = 0.0;
    }
 
    if (diag != nullptr) {
        diag->P_elastic = P_el;
        diag->P_inelastic = P_inel;
        diag->P_chemical = P_chem;
        diag->P_chemical_e = P_chem_e;
        diag->P_Joule = P_joule;
        // Collision-frequency diagnostics (Mitchner II-13.3 / II-8.11e /
        // VIII-3.8), from the same transport() call that produced sigma/nu_E
        // for this state.
        diag->nu_eH = nu_eH;
        diag->nu_ee = CollisionModel::electronElectronCollisionFrequency(n_e, Te_c);
        diag->nu_eI = nu_eI;
        diag->nu_eH_mass_weighted = nu_eH_mass_weighted;
    }
}
 
double ReactorRHS::conductivity(double Tg, double Te, const double* Y, double rho) const
{
    double P, Xe, Tmean, n_e;
    setState(Tg, Te, Y, rho, P, Xe, Tmean, n_e);
    double sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted;
    transport(std::max(Te, 200.0), n_e, P, Tmean, sigma, nu_E, nu_eH, nu_eI,
              nu_eH_mass_weighted);
    return sigma;
}
 
double ReactorRHS::electronDensity(double Tg, double Te, const double* Y,
                                    double rho) const
{
    double P, Xe, Tmean, n_e;
    setState(Tg, Te, Y, rho, P, Xe, Tmean, n_e);
    return n_e;
}
 
} // namespace rizer