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.plasma.plasma_extension._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_dU.resize(nR);
    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);
    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;
        }
    }
}
 
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& cond_maxwell_2) 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, cond_maxwell_2);
    } 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). cond_maxwell_2 needs 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;
        cond_maxwell_2 = 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, cond_maxwell_2;
    transport(Te_c, n_e, P, Tmean, sigma, nu_E, cond_maxwell_2);
 
    // 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);
        // Internal-energy change per reaction (constant-volume reactor), matching
        // the Python reference: @f$dU_i = \sum_k \nu_{ki} u_k@f$ =
        // getReactionDelta(u). Using dU (not dH) keeps the R*T_k pressure-work of
        // created particles out of the electron->gas inelastic exchange (see
        // plasma_power_inelastic).
        kin->getReactionDelta(Cantera::span<const double>(m_uk.data(), m_nsp),
                              Cantera::span<double>(m_dU.data(), m_dU.size()));
        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;
        for (std::size_t i = 0; i < m_rop.size(); i++) {
            if (m_is_electron_rxn[i]) P_inel += m_dU[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;
        // Maxwellian-distribution conditions (Mitchner VIII-3.8): C1 = nu_ee/m_e,
        // C2 = sum_h nu_eh/m_h (cond_maxwell_2 above, from the same transport()
        // call that produced sigma/nu_E for this state).
        const double nu_ee = CollisionModel::electronElectronCollisionFrequency(n_e, Te_c);
        diag->cond_maxwell_1 = nu_ee / units::m_e;
        diag->cond_maxwell_2 = cond_maxwell_2;
    }
}
 
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, cond_maxwell_2;
    transport(std::max(Te, 200.0), n_e, P, Tmean, sigma, nu_E, cond_maxwell_2);
    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