PlasmaReactor0D.cpp

Go to the documentation of this file.

#include "PlasmaReactor0D.h"
#include "numeric/Format.h"
#include "kinetics/PlasmaRates.h"
#include "thermo/PlasmaState.h"
#include "units.h"
 
#include "cantera/thermo/PlasmaPhase.h"
#include "cantera/thermo/ThermoPhase.h"
 
#include <algorithm>
#include <cmath>
#include <sstream>
#include <stdexcept>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   Constructor                                                                  //
    // ------------------------------------------------------------------------------ //
 
    PlasmaReactor0D::PlasmaReactor0D(const PlasmaReactor0DConfig &cfg)
        : m_sol(loadMechanism(cfg.mech, cfg.phase)), m_gap(cfg.gap)
    {
        auto thermo = m_sol->thermo();
        const std::size_t nsp = thermo->nSpecies();
        m_plasma = dynamic_cast<Cantera::PlasmaPhase *>(thermo.get());
        if (m_plasma == nullptr)
        {
            throw std::invalid_argument(
                "PlasmaReactor0D: phase \"" + cfg.phase + "\" in \"" + cfg.mech +
                "\" is a " + thermo->type() +
                ", not a PlasmaPhase, so it carries no electron temperature.");
        }
        m_ie = m_plasma->electronSpeciesIndex();
        // 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
        // Refused rather than resized: padding gives the missing species no
        // collision model and truncating drops the caller's tail, so either
        // way a species silently leaves every collision sum.
        if (cfg.specs.size() != nsp)
        {
            throw std::invalid_argument(
                "PlasmaReactor0D: got collision specs for " +
                std::to_string(cfg.specs.size()) + " species, but \"" + cfg.mech +
                "\" has " + std::to_string(nsp) +
                " species. Give one spec per mechanism species, in its order.");
        }
        std::vector<CollisionModel::SpeciesSpec> specs = cfg.specs;
        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.
        CollisionModel collision(specs, cfg.Te_min, cfg.Te_max, cfg.Te_n);
 
        // PlasmaCell is the physics core shared with PlasmaChannel1D: chemistry source
        // terms (if reacting), Joule heating, and elastic/inelastic Tg<->Te exchange.
        PlasmaCellConfig rc;
        rc.reacting = cfg.reacting;
        rc.collision = std::move(collision);
        rc.species_h0k = cfg.species_h0k;
        m_rhs = std::make_unique<PlasmaCell>(m_sol, std::move(rc));
    }
 
    // ------------------------------------------------------------------------------ //
    //   Public methods                                                               //
    // ------------------------------------------------------------------------------ //
 
    std::vector<std::string> PlasmaReactor0D::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 PlasmaReactor0D::rhs(double Tg, double Te, double V, const double *Y, double E,
                              double mass, double p_ext, double polytropic_index,
                              double *dydt) const
    {
        // Three steps:
        //   1. the shared cell's source terms, at the density the volume implies
        //   2. the polytropic volume equation, and the pressure work it does on
        //      each pool -- skipped entirely at constant volume
        //   3. pack dydt = [dTg, dTe, dV, dY_0..dY_{K-1}]
 
        // ---- 1. the shared cell at this state -------------------------------
        // 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$.
        // +inf is the constant-volume case. Any other index has to be a usable
        // polytropic one: running it at constant volume would drop the
        // requested expansion without a word. Checked before the source
        // terms, so a setup error is not hidden behind a recoverable
        // TrialStateError from a refused state.
        const bool constant_volume = std::isinf(polytropic_index) && polytropic_index > 0.0;
        if (!constant_volume && !(polytropic_index > 0.0))
        {
            throw std::invalid_argument(
                "PlasmaReactor0D::rhs: polytropic_index must be positive, or +inf for a "
                "constant volume; got " + significant(polytropic_index) + ".");
        }
        if (!constant_volume && !(m_gap > 0.0))
        {
            throw std::invalid_argument(
                "PlasmaReactor0D::rhs: a finite polytropic_index needs gap > 0, since "
                "the rarefaction timescale comes from that length; got gap = " +
                significant(m_gap) + " m.");
        }
        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->sourceTerms(Tg, Te, Y, rho, E, dTg, dTe, dydt + 3, lhsTg, lhsTe, &m_diag);
 
        // ---- 2. polytropic volume, and its pressure work --------------------
        double dV = 0.0;
        if (!constant_volume)
        {
            // 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 TwoTemperatureState eos = readTwoTemperatureState(*m_plasma);
            const double P = eos.pressure;
            // Electron partial pressure [Pa], PHYSICS.md (G6): @f$p_e = n_e k_B T_e@f$,
            // taken from PlasmaPhase rather than as the Dalton share @f$X_e P@f$. The
            // two differ by @f$\overline{T}/T_e@f$, because Cantera's P is
            // @f$n_\text{tot}k_B\overline{T}@f$ -- a factor ~10 at
            // @f$T_e = 10 T_g@f$, all of it in the pressure work below.
            const double p_e = eos.electron_pressure;
            // The temperatures need no check of their own here: step 1's
            // sourceTerms() already refused, as a recoverable TrialStateError,
            // any Tg outside [TemperatureFloor, TemperatureCeiling] and any Te
            // under the floor. Te above the ceiling is valid: the electron's
            // cp/R is the constant 2.5, so the mixture cp and cv below hold.
            // Two-temperature mixture temperature, the effective temperature of
            // the speed-of-sound estimate below. Read from the phase rather than
            // rebuilt: PlasmaPhase::meanTemperature() is
            // @f$\overline{T} = (1 - X_e) T_g + X_e T_e@f$ by definition, and it
            // is the same mean the `pressure` above was computed with.
            const double Tmean = eos.mean_temperature;
            const double cp = thermo->cp_mass();
            const double cv = thermo->cv_mass();
            // A gas mixture has c_v > 0; anything else means the state handed to
            // the thermo layer is not one, and the speed of sound below has no
            // meaning. Raised rather than substituted: a stand-in gamma would
            // return a plausible rarefaction rate computed from nothing.
            if (!(cv > 0.0))
            {
                std::ostringstream message;
                message << "PlasmaReactor0D::rhs: heat capacity at constant "
                           "volume is "
                        << std::scientific << cv
                        << " J/kg/K at Tg=" << Tg << " K, Te=" << Te
                        << " K, rho=" << rho
                        << " kg/m^3. It must be positive for the polytropic "
                           "volume equation's speed of sound. Check the "
                           "composition and the temperatures at this state.";
                throw std::runtime_error(message.str());
            }
            const double gamma = cp / cv;
            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;
            // Raised, never zeroed, for the reason PlasmaCell::sourceTerms gives:
            // a substituted 0 hands the integrator a smooth derivative it accepts,
            // and the volume would freeze on a state the solver never asked about.
            if (!std::isfinite(dV))
            {
                std::ostringstream message;
                message << "PlasmaReactor0D::rhs: non-finite volume rate dV/dt = "
                        << dV << " m^3/s at P = " << P << " Pa, p_ext = " << p_ext
                        << " Pa, V = " << V << " m^3, tau = " << tau
                        << " s, polytropic index " << polytropic_index << ".";
                throw std::runtime_error(message.str());
            }
            // Pressure work removes energy from each population as the cell expands.
            dTe -= (p_e * dV / V) / lhsTe;
            dTg -= ((P - p_e) * dV / V) / lhsTg;
        }
 
        // ---- 3. pack the state derivative -----------------------------------
        dydt[0] = dTg;
        dydt[1] = dTe;
        dydt[2] = dV;
        // dydt[3..] (dY) already written by rates().
    }
 
} // namespace rizer