PlasmaCell.cpp

Go to the documentation of this file.

#include "PlasmaCell.h"
#include "numeric/TableLookup.h"
 
#include "kinetics/ElectronReactions.h"
#include "thermo/PlasmaState.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 <stdexcept>
#include <string>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   Constructor                                                                  //
    // ------------------------------------------------------------------------------ //
 
    PlasmaCell::PlasmaCell(std::shared_ptr<Cantera::Solution> sol,
                           PlasmaCellConfig cfg)
        : m_sol(std::move(sol)),
          m_nsp(m_sol->thermo()->nSpecies()),
          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, "sigma_T/sigma_v")),
          m_nu_m_tab(PropertyTable::fromTableOrScalar(
              cfg.nu_m_T, cfg.nu_m_v, "nu_m_T/nu_m_v")),
          m_nu_m(cfg.nu_m),
          m_nu_E(cfg.nu_E)
    {
        m_plasma = dynamic_cast<Cantera::PlasmaPhase *>(m_sol->thermo().get());
        if (m_plasma == nullptr)
        {
            throw std::invalid_argument(
                "PlasmaCell: the phase is a " + m_sol->thermo()->type() +
                ", not a PlasmaPhase, so it carries no electron temperature.");
        }
        m_ie = m_plasma->electronSpeciesIndex();
        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);
        // Which reactions the electron takes part in and what each costs it:
        // a property of the MECHANISM, built in kinetics/ElectronReactions.h.
        const ElectronReactionIndex electron_reactions =
            buildElectronReactionIndex(*m_sol->kinetics(),
                                       m_sol->thermo()->speciesName(m_ie),
                                       cfg.species_h0k, m_nsp, m_reacting);
        m_is_electron_rxn = std::move(electron_reactions.is_electron_reaction);
        m_eps_th = std::move(electron_reactions.threshold_energy);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Public methods                                                               //
    // ------------------------------------------------------------------------------ //
 
    std::string PlasmaCell::speciesName(std::size_t k) const
    {
        return m_sol->thermo()->speciesName(k);
    }
 
    void PlasmaCell::sourceTerms(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
    {
        // Seven steps, in order:
        //   1. reject a state outside the range the fits hold over
        //   2. put the state on the phase, read the two-temperature EOS back
        //   3. capacitances lhs_Tg, lhs_Te -- the J/m^3/K the balances divide by
        //   4. transport closure: sigma and the collision frequencies
        //   5. the field-driven powers, Joule and elastic exchange
        //   6. chemistry: dY, the two chemical powers and the inelastic one
        //   7. assemble both balances, check them, fill the diagnostics
 
        auto thermo = m_sol->thermo();
        auto kin = m_sol->kinetics();
 
        // ---- 1-2. reject an invalid state, set it, read the EOS back --------
        // setState() raises a TrialStateError for a temperature under
        // TemperatureFloor, so CVODE rejects the trial step and retries a
        // smaller one rather than abandoning the solve.
        const TwoTemperatureState state = setState(Tg, Te, Y, rho);
        const double n_e = state.electron_density;
 
        // ---- 3. energy-equation capacitances --------------------------------
        // Also leaves m_cpR current, which the heavy-pool capacitance below reads.
        fillSpeciesCp();
        const double We = m_wt[m_ie];
        const double cve = monatomicElectronCv(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;
        // Raised, never floored: the gas balance divides by it, and a stand-in
        // returns a finite rate for a state with no heavy gas to heat. With
        // the temperatures inside the fits every heavy cv is positive, so this
        // fires only on a composition with no heavy species in it.
        lhs_Tg = cv_heavy * rho;
        if (!(lhs_Tg > 0.0))
        {
            throw TrialStateError(
                "PlasmaCell::sourceTerms",
                "the heavy-pool heat capacity is " + std::to_string(lhs_Tg) +
                " J/m^3/K at Tg = " + std::to_string(Tg) + " K, Te = " +
                std::to_string(Te) + " K, rho = " + std::to_string(rho) +
                " kg/m^3; it must be positive for the gas-temperature balance.");
        }
 
        // ---- 4. transport closure -------------------------------------------
        double sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted;
        transport(Te, state, sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted);
 
        // ---- 5. field-driven powers -----------------------------------------
        // 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 - Tg);
        // 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;
 
        // ---- 6. chemistry: dY and the chemical powers -----------------------
        double P_chem_H = 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_H += m_uk[k] * m_wdot[k];
            }
            P_chem_e = m_wdot[m_ie] * cve * We * Te;
            // 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];
            }
        }
 
        // ---- 7. assemble the two balances -----------------------------------
        dTe = (-P_chem_e - P_el - P_inel + P_joule) / lhs_Te;
        dTg = (-P_chem_H + P_el + P_inel) / lhs_Tg;
 
        // Raised, never zeroed. A non-finite rate means the state this was
        // evaluated at is not one the closures can describe -- an overflowing
        // rate constant, or a trial composition the integrator should back away
        // from. Substituting 0.0 hands CVODE a smooth, well-behaved derivative
        // instead, so its error test passes and the step is ACCEPTED: the
        // integration then continues from the corrupted state and returns a
        // plausible wrong answer rather than failing. A TrialStateError lets
        // CVODE reject the step and retry smaller, which is what its error
        // control is for.
        if (!std::isfinite(dTg) || !std::isfinite(dTe))
        {
            throw TrialStateError(
                "PlasmaCell::sourceTerms",
                "non-finite temperature rate at Tg = " +
                std::to_string(Tg) +
                " K, Te = " +
                std::to_string(Te) +
                " K, rho = " +
                std::to_string(rho) +
                " kg/m^3 (dTg = " +
                std::to_string(dTg) +
                ", dTe = " +
                std::to_string(dTe) + ").");
        }
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            if (!std::isfinite(dY[k]))
            {
                throw TrialStateError(
                    "PlasmaCell::sourceTerms",
                    "non-finite production rate for species " +
                    std::to_string(k) +
                    " at Tg = " +
                    std::to_string(Tg) +
                    " K, Te = " +
                    std::to_string(Te) +
                    " K, rho = " +
                    std::to_string(rho) +
                    " kg/m^3.");
            }
        }
 
        // Fill the optional Diagnostics struct if the caller passed a non-null pointer.
        if (diag != nullptr)
        {
            diag->P_elastic = P_el;
            diag->P_inelastic = P_inel;
            diag->P_chemical_H = P_chem_H;
            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);
            diag->nu_eI = nu_eI;
            diag->nu_eH_mass_weighted = nu_eH_mass_weighted;
        }
    }
 
    double PlasmaCell::electricalConductivity(double Tg, double Te, const double *Y,
                                              double rho) const
    {
        const TwoTemperatureState state = setState(Tg, Te, Y, rho);
        double sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted;
        transport(Te, state, sigma, nu_E, nu_eH, nu_eI, nu_eH_mass_weighted);
        return sigma;
    }
 
    double PlasmaCell::electronDensity(double Tg, double Te, const double *Y,
                                       double rho) const
    {
        return setState(Tg, Te, Y, rho).electron_density;
    }
 
    void PlasmaCell::pressureAndElectronDensity(double Tg, double Te,
                                                const double *Y, double rho,
                                                double &p,
                                                double &electron_density) const
    {
        const TwoTemperatureState state = setState(Tg, Te, Y, rho);
        p = state.pressure;
        electron_density = state.electron_density;
    }
 
    void PlasmaCell::updateSpeciesCp(double Tg, double Te, const double *Y,
                                     double rho) const
    {
        setState(Tg, Te, Y, rho);
        fillSpeciesCp();
    }
 
    void PlasmaCell::speciesEnthalpies(double Tg, double Te, const double *Y,
                                       double rho, double *h) const
    {
        setState(Tg, Te, Y, rho);
        m_sol->thermo()->getPartialMolarEnthalpies(
            Cantera::span<double>(h, m_nsp));
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            h[k] /= m_wt[k];
        }
        const double We = m_wt[m_ie];
        h[m_ie] = (monatomicElectronCv(We) + GasConstant / We) * Te;
    }
 
    // ------------------------------------------------------------------------------ //
    //   Private methods                                                              //
    // ------------------------------------------------------------------------------ //
 
    void PlasmaCell::fillSpeciesCp() const
    {
        // 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, for consistency with
        // how cve/lhs_Te bypass Cantera's own electron Cp in sourceTerms().
        m_sol->thermo()->getCp_R(m_cpR);
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            m_cp_mass[k] = m_cpR[k] * GasConstant / m_wt[k];
        }
        const double We = m_wt[m_ie];
        m_cp_mass[m_ie] = monatomicElectronCv(We) + GasConstant / We;
    }
 
    TwoTemperatureState PlasmaCell::setState(double Tg, double Te, const double *Y,
                                             double rho) const
    {
        // Rejected rather than clamped. Every rate, capacity and collision
        // frequency is a fit that stops at TemperatureFloor, so a colder state
        // is a diverged iterate, and clamping it would hand the caller a value
        // for a temperature it never proposed. Tg is also held under
        // TemperatureCeiling, above which the heavy-species thermo fits
        // extrapolate: the C+, C++ and C+++ fits of the CH4 mechanism give
        // cp/R < 1 from 1.52e5 K, i.e. a negative heavy-pool heat capacity.
        // Te has no such ceiling here: the electron's own heat capacity is the
        // analytic monatomic value, and the P4 and pulsed-reactor runs reach
        // Te above 1e5 K. Checked here, where every public read sets its
        // state, so no read can skip it. The negated comparisons also reject a
        // NaN, which `Tg < TemperatureFloor` would wave through.
        if (!(Tg >= TemperatureFloor) || !(Te >= TemperatureFloor) ||
            !(Tg <= TemperatureCeiling))
        {
            std::ostringstream message;
            message << "Tg = " << Tg << " K, Te = "
                    << Te << " K, rho = " << rho
                    << " kg/m^3; Tg must lie within [" << TemperatureFloor << ", "
                    << TemperatureCeiling << "] K and Te must be at least "
                    << TemperatureFloor
                    << " K, where the mechanism's thermo and rate fits hold.";
            throw TrialStateError("PlasmaCell::setState", message.str());
        }
        auto thermo = m_sol->thermo();
        // Negative mass fractions are clipped to zero before Cantera sees them;
        // setMassFractions_NoNorm keeps them un-renormalized so Y still sums to
        // (numerically) 1 as the integrator expects. Clipped rather than passed
        // through, as Cantera's own reactors and the Python reference do, because
        // the closures cannot take a negative density: a negative n_e has no
        // Coulomb logarithm. Measured over the suite and the native examples: the
        // native channel's negatives stay at round-off (1e-13 or below, trial
        // states), held there by CVODE's tolerances, since Cantera's CVODE wrapper
        // applies no positivity constraint; the Python split solver's scipy vode
        // has none either and reaches -1.3e-7 (C2H6+). A NaN is
        // a broken state rather than round-off, so it raises instead of being
        // clipped to zero.
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            if (std::isnan(Y[k]))
            {
                throw TrialStateError(
                    "PlasmaCell::setState",
                    "the mass fraction of species " + m_sol->thermo()->speciesName(k) +
                    " is NaN at Tg = " + std::to_string(Tg) + " K, Te = " +
                    std::to_string(Te) + " K.");
            }
            m_Yc[k] = (Y[k] > 0.0) ? Y[k] : 0.0;
        }
        thermo->setMassFractions_NoNorm(Cantera::span<const double>(m_Yc.data(), m_nsp));
        thermo->setElectronTemperature(Te);
        thermo->setState_TD(Tg, rho);
        // Read back rather than rebuilt: the mean temperature and the electron
        // density follow from the state just set, and PlasmaPhase is the one
        // that defines them (thermo/PlasmaState.h). Returned whole rather than
        // unpacked into out-parameters, so a caller takes the one field it
        // wants instead of declaring four locals to drop three of them.
        return readTwoTemperatureState(*m_plasma);
    }
 
    void PlasmaCell::transport(double Te, const TwoTemperatureState &state,
                               double &sigma, double &nu_E, double &nu_eH,
                               double &nu_eI, double &nu_eH_mass_weighted) const
    {
        const double n_e = state.electron_density;
        // 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 PlasmaCellConfig in PlasmaCell.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 = state.pressure / (Boltzmann * state.mean_temperature);
            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, 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
                    : evalAtState(m_nu_m_tab, Te, "PlasmaCell::transport", "nu_m(Te)", "K",
                                  "Extend the nu_m_T/nu_m_v table.");
            if (nu_m_eff > 0.0)
            {
                sigma = n_e * units::e * units::e / (units::m_e * nu_m_eff);
            }
            else
            {
                sigma = m_sigma.empty()
                            ? 0.0
                            : evalAtState(m_sigma, Te, "PlasmaCell::transport",
                                          "sigma(Te)", "K",
                                          "Extend the sigma_T/sigma_v table.");
            }
            nu_E = m_nu_E;
            nu_eH = 0.0;
            nu_eI = 0.0;
            nu_eH_mass_weighted = 0.0;
        }
    }
 
} // namespace rizer