PlasmaChannel1D.h

Go to the documentation of this file.

 
#ifndef RIZER_PLASMA_CHANNEL_1D_H
#define RIZER_PLASMA_CHANNEL_1D_H
 
#ifndef NOMINMAX
#define NOMINMAX
#endif
#include "cantera/oneD/Domain1D.h"
#include "cantera/base/Solution.h"
#include "transport/CollisionModel.h"
#include "numeric/PropertyTable.h"
#include "BoundaryConditions.h"
#include "PlasmaCell.h"
#include "RadialDiffusion.h"
 
#include <memory>
#include <string>
#include <vector>
 
namespace rizer
{
 
    // Physics + initial/boundary configuration for a PlasmaChannel1D.
    struct PlasmaChannelConfig
    {
        // ---- Mechanism ----------------------------------------------------

        std::string mech;
        std::string phase;
 
        // ---- Grid ---------------------------------------------------------

        std::size_t n_points = 1;
        double R_max = 1.0e-3;
        std::vector<double> grid;
 
        // ---- Initial state ------------------------------------------------

        double rho = 0.1;
        std::vector<double> initTg_r, initTg_v;
        std::vector<double> initTe_r, initTe_v;
        std::vector<double> initY_r, initY_v;
        std::vector<double> init_rho_r, init_rho_v;
        std::vector<double> init_ur_r, init_ur_v;
 
        // ---- Chemistry ----------------------------------------------------

        bool reacting = true;
        std::vector<double> species_h0k;
 
        // ---- Discharge ----------------------------------------------------

        double electric_field = 0.0;
 
        // ---- Transport, bulk closures -------------------------------------
        //
        // Each `_T` array is a table's temperature axis [K] -- electron
        // temperature for sigma and kappa_e, heavy-pool temperature for
        // kappa_H -- and each `_v` array the tabulated value there.

        std::vector<double> sigma_T, sigma_v;
        double nu_E = 0.0;
        double nu_m = 0.0;
        std::vector<double> nu_m_T, nu_m_v;
        std::vector<double> kappa_H_T, kappa_H_v;
        std::vector<double> kappa_e_T, kappa_e_v;
        double D_species = 0.0;
 
        // ---- Transport, composition-resolved model ------------------------

        std::vector<CollisionModel::SpeciesSpec> specs;
        double Te_min = 200.0;
        double Te_max = 1.0e5;
        std::size_t Te_n = 400;
 
        // ---- Ambient and boundary -----------------------------------------

        double T_amb = 300.0;
        double p_amb = -1.0;
        std::vector<double> Y_amb;
        std::string outer_bc_mode = "ambient_ghost";
 
        // ---- Hydrodynamics ------------------------------------------------

        bool hydrodynamics = false;
        double cfl = 0.9;
        double hydro_substep = 0.0;
    };
 
    class PlasmaChannel1D : public Cantera::Domain1D
    {
    public:
        explicit PlasmaChannel1D(const PlasmaChannelConfig &cfg);
 
        std::string domainType() const override { return "plasma-channel-1d"; }

        void eval(std::size_t jg, Cantera::span<const double> xg, Cantera::span<double> rg,
                  Cantera::span<int> maskg, double rdt) override;

        std::string componentName(std::size_t n) const override;

        std::size_t componentIndex(const std::string &name,
                                   bool checkAlias = true) const override;
        bool hasComponent(const std::string &name, bool checkAlias = true) const override;

        void getValues(const std::string &component,
                       Cantera::span<double> values) const override;

        double initialValue(std::size_t n, std::size_t j) override;

        void setElectricField(double electricField) { m_electric_field = electricField; }
        std::size_t nSpecies() const { return m_n_species; }

        const PropertyTable &heavyConductivity() const { return m_kappa_H; }

        const PropertyTable &electronConductivity() const { return m_kappa_e; }

        const std::vector<double> &molecularWeights() const
        {
            return m_cell->molecularWeights();
        }

        const std::vector<double> &speciesCp() const { return m_cell->speciesCp(); }

        void updateSpeciesCp(double Tg, double Te, const double *Y, double rho) const
        {
            m_cell->updateSpeciesCp(Tg, Te, Y, rho);
        }

        OuterBc outerBcMode() const { return m_outer_bc_mode; }

        bool hasTransport() const { return m_has_transport; }

        bool hydrodynamics() const { return m_hydrodynamics; }

        std::size_t electronIndex() const { return m_electron_index; }

        const std::vector<double> &density() const { return m_rho; }

        const std::vector<double> &velocity() const { return m_u_r; }

        void setHydroState(const std::vector<double> &rho,
                           const std::vector<double> &u_r);

        void nodeDiagnostics(double Tg, double Te, const double *Y, std::size_t j,
                             double &pressure, double &electron_density) const;

        void nodeRates(std::size_t j, double Tg, double Te, const double *Y,
                       double electric_field, double *ydot) const;

        double nodeElectricalConductivity(double Tg, double Te, const double *Y,
                                          std::size_t j) const;

        void wallEnergyOutflow(Cantera::span<const double> x, double &conduction_out,
                               double &diffusion_out,
                               std::vector<double> &species_mass_out) const;

        const std::vector<double> &ambientMassFractions() const { return m_Y_amb; }

        double faceRadius(std::size_t j) const;

        double faceSpacing(std::size_t j) const;

        void setEvaluatorPoolSize(std::size_t slots);

        std::size_t evaluatorPoolSize() const { return 1 + m_extra_rhs.size(); }
 
    private:
        RadialDiffusionCoefficients diffusionCoefficients() const;
 
        // Cantera solution, which owns the thermo/kinetics manager (PlasmaPhase) and the
        // transport manager (Transport).
        std::shared_ptr<Cantera::Solution> m_solution;
        // The single evaluator slot 0, which owns its own Solution and is the one
        // nodeRates() uses when the pool size is 1. It is also the one that eval()
        // uses, so it is the only one that sees the evolving state vector.
        std::unique_ptr<PlasmaCell> m_cell;
 
        // Evaluator slots 1..n for setEvaluatorPoolSize(); slot 0 is
        // m_cell. Each owns its own Solution, so nothing is shared.
        std::vector<std::unique_ptr<PlasmaCell>> m_extra_rhs;
        // What setEvaluatorPoolSize() needs to build another evaluator: the same
        // closure configuration, over a mechanism loaded again from scratch.
        PlasmaCellConfig m_reactor_config;
        std::string m_mech, m_phase;
        std::size_t m_n_species;            // Number of species K [-] (thermo->nSpecies())
        std::size_t m_electron_index;       // Species index of "e-" [-] (npos if the
                                            // mechanism has none)
        double m_outer_radius;              // Outer (wall) radius R [m]
        std::vector<double> m_rho;          // Mass density rho [kg/m^3], one per radial
                                            // node. Held outside the Domain1D solution
                                            // vector, which stays nv=2+K: the radial
                                            // hydrodynamic sub-step owns this field
                                            // (NUMERICS.md section 4.2), so eval() has
                                            // no residual row to write for it
        std::vector<double> m_u_r;          // Radial velocity u_r [m/s], one per radial
                                            // node, owned by the hydro sub-step like
                                            // m_rho and held outside the solution
                                            // vector for the same reason
        bool m_hydrodynamics;               // Run the radial hydrodynamic sub-step [-]
        OuterBc m_outer_bc_mode;            // Outer boundary treatment for the hydro
                                            // sub-step, parsed from the config string
        double m_wall_ghost_density;        // Mass density of the wall ghost half-cell
                                            // [kg/m^3]. The ghost carries the ambient
                                            // composition and temperature, so with
                                            // hydrodynamics on -- the only mode that
                                            // supplies an ambient pressure -- it carries
                                            // the ambient density too, rather than the
                                            // outermost node's
        double m_electric_field;            // Applied electric field E [V/m], set per
                                            // eval() via setElectricField()
        double m_species_diffusivity;       // Constant scalar species mass diffusivity
                                            // D [m^2/s] (0 = frozen)
        bool m_reacting;                    // Include finite-rate chemistry source terms [-]
        double m_ambient_temperature;       // Wall/ambient temperature T_amb [K]
        std::vector<double> m_Y_amb;        // Ambient mass fractions [-] (length K,
                                            // sums to 1), outside the domain
        PropertyTable m_kappa_H;            // Heavy-pool thermal conductivity table
                                            // @f$\kappa_\mathrm{H}(T_g)@f$ [W/m/K];
                                            // empty -> disabled
        PropertyTable m_kappa_e;            // Electron thermal conductivity table
                                            // @f$\kappa_e(T_e)@f$ [W/m/K]; empty ->
                                            // disabled
        PropertyTable m_initTg;             // Initial radial profile Tg(r) [K]; empty
                                            // -> uniform m_ambient_temperature
        PropertyTable m_initTe;             // Initial radial profile Te(r) [K]; empty
                                            // -> uniform m_ambient_temperature
        std::vector<PropertyTable> m_initY; // Initial radial profiles Y_k(r)
                                            // [-], one per species
        bool m_has_transport;               // True if any conduction/diffusion term is
                                            // enabled (kappa_H, kappa_e, or
                                            // D_species > 0)
    };
 
} // namespace rizer
 
#endif