Plasma0DReactor.cpp

Go to the documentation of this file.

#include "Plasma0DReactor.h"
#include "kinetics/PlasmaRates.h"
#include "units.h"
 
#include "cantera/thermo/ThermoPhase.h"
 
#include <algorithm>
#include <cmath>
 
namespace rizer {
 
namespace {
Cantera::shared_ptr<Cantera::Solution>
loadMechanism(const std::string& mech, const std::string& phase)
{
    registerPlasmaRates();   // Idempotent
    return Cantera::newSolution(mech, phase);
}
} // namespace
 
Plasma0DReactor::Plasma0DReactor(const Config& cfg)
    : m_sol(loadMechanism(cfg.mech, cfg.phase))
    , m_gap(cfg.gap)
{
    auto thermo = m_sol->thermo();
    const std::size_t nsp = thermo->nSpecies();
    // Locate the electron pseudo-species; its index is used throughout (mole/mass
    // fraction lookups, is_electron flag below, and the partial-pressure split in
    // rhs()).
    m_ie = thermo->speciesIndex("e-");
    if (m_ie == Cantera::npos) {
        throw std::invalid_argument("Plasma0DReactor: mechanism has no 'e-' species.");
    }
    // Fill molar_mass [kg/mol] and is_electron into the per-species collision specs.
    std::vector<double> wt(thermo->molecularWeights().begin(),
                           thermo->molecularWeights().end());   // kg/kmol
    std::vector<CollisionModel::SpeciesSpec> specs = cfg.specs;
    specs.resize(nsp);
    for (std::size_t k = 0; k < nsp; k++) {
        specs[k].molar_mass = wt[k] * 1.0e-3;   // kg/kmol -> kg/mol
        specs[k].is_electron = (k == m_ie);
    }
    // CollisionModel tabulates composition-resolved transport/rate coefficients
    // (e.g. mobility, collision frequency) over the electron-temperature grid
    // [Te_min, Te_max] with Te_n points; spitzer toggles the Coulomb contribution.
    CollisionModel collision(specs, cfg.Te_min, cfg.Te_max, cfg.Te_n, cfg.spitzer);
 
    // ReactorRHS is the physics core shared with PlasmaChannel1D: chemistry source
    // terms (if reacting), Joule heating, and elastic/inelastic Tg<->Te exchange.
    ReactorRHS::Config rc;
    rc.reacting = cfg.reacting;
    rc.collision = std::move(collision);
    rc.species_h0k = cfg.species_h0k;
    m_rhs = std::make_unique<ReactorRHS>(m_sol, std::move(rc));
}
 
std::vector<std::string> Plasma0DReactor::speciesNames() const
{
    std::vector<std::string> names(m_rhs->nSpecies());
    for (std::size_t k = 0; k < names.size(); k++) names[k] = m_rhs->speciesName(k);
    return names;
}
 
void Plasma0DReactor::rhs(double Tg, double Te, double V, const double* Y, double E,
                           double mass, double p_ext, double polytropic_index,
                           double* dydt) const
{
    // Constant-mass reactor: density is derived from the instantaneous volume rather
    // than integrated directly (mirrors Isomass2TVolumeReactor in Python).
    const double rho = mass / V;
    // dTg, dTe accumulate the RHS's energy-balance source terms [K/s]; lhsTg,
    // lhsTe are the (heat-capacity-like) coefficients multiplying dTg/dt, dTe/dt
    // on the LHS of the energy balance, so source terms below can be added as
    // @f$ \text{source} / \text{lhs} @f$.
    double dTg = 0.0, dTe = 0.0, lhsTg = 1.0, lhsTe = 1.0;
    // dY goes straight into dydt[3..]; rates() leaves thermo at (Tg,Te,Y,rho).
    // &m_diag caches the power/Maxwellian-condition breakdown for diagnostics().
    m_rhs->rates(Tg, Te, Y, rho, E, dTg, dTe, dydt + 3, lhsTg, lhsTe, &m_diag);
 
    double dV = 0.0;
    if (std::isfinite(polytropic_index) && polytropic_index > 0.0 && m_gap > 0.0) {
        // Polytropic rarefaction (mirrors Isomass2TVolumeReactor): the hot
        // core's overpressure relaxes toward p_ext over the acoustic (rarefaction)
        // timescale tau derived below.
        auto thermo = m_rhs->solution()->thermo();
        const double P = thermo->pressure();
        const double Xe = thermo->moleFraction(m_ie);
        // Electron partial pressure [Pa], via the ideal-gas Dalton split
        // @f$ p_e = X_e P @f$.
        const double p_e = Xe * P;
        // Clamp to avoid unphysical/singular values (e.g. at startup) feeding
        // sqrt() below.
        const double Tg_c = std::max(Tg, 200.0);
        const double Te_c = std::max(Te, 200.0);
        // Two-temperature mixture temperature (mole-fraction-weighted blend of Tg
        // and Te), used as the effective temperature for the speed-of-sound
        // estimate below:
        // @f[
        //     \overline{T} = (1 - X_e) T_g + X_e T_e.
        // @f]
        const double Tmean = (1.0 - Xe) * Tg_c + Xe * Te_c;
        const double cp = thermo->cp_mass();
        const double cv = thermo->cv_mass();
        // Fall back to diatomic-air @f$ \gamma @f$ if @f$ c_v \approx 0 @f$.
        const double gamma = (cv > 0.0) ? cp / cv : 1.4;
        const double M = thermo->meanMolecularWeight() * 1.0e-3;   // kg/mol
        // Treat the plasma volume as a cylinder of fixed axial length m_gap to
        // recover an effective radius r, then the ideal-gas speed of sound at
        // @f$ \overline{T} @f$ and the acoustic (rarefaction) timescale tau over
        // which the core relaxes toward p_ext:
        // @f[
        //     r = \sqrt{\frac{V}{\text{m\_gap} \, \pi}}, \qquad
        //     c_\text{sound} = \sqrt{\frac{\gamma R \overline{T}}{M}}, \qquad
        //     \tau = \frac{r}{c_\text{sound}}.
        // @f]
        const double radius = std::sqrt(V / (m_gap * units::pi));
        const double c_sound = std::sqrt(gamma * units::R * Tmean / M);
        const double tau = radius / c_sound;
        // Polytropic volume equation (mirrors Isomass2TVolumeReactor): the
        // pressure relaxes toward p_ext over tau, and the volume responds via the
        // polytropic relation with index n:
        // @f[
        //     \frac{dp}{dt} = -\frac{P - p_\text{ext}}{\tau}, \qquad
        //     \frac{dV}{dt} = -\frac{1}{n} \frac{V}{P} \frac{dp}{dt}.
        // @f]
        const double dp_dt = -(P - p_ext) / tau;
        dV = -1.0 / polytropic_index * (V / P) * dp_dt;
        // Pressure work removes energy from each population as the cell expands.
        if (std::isfinite(dV)) {
            dTe -= (p_e * dV / V) / lhsTe;
            dTg -= ((P - p_e) * dV / V) / lhsTg;
        } else {
            dV = 0.0;
        }
    }
 
    dydt[0] = dTg;
    dydt[1] = dTe;
    dydt[2] = dV;
    // dydt[3..] (dY) already written by rates().
}
 
} // namespace rizer