diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index 1cddd48df41..7b8508f6c00 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -211,6 +211,18 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! Update the value of #m_crossSectionsInterpolated [m2] void updateInterpolatedCrossSection(span); + //! Returns the name of the collision linking the reaction to the data stored in electron-collision + const string& collisionName() const; + + //! Return whether this rate object has tabulated cross-section data. + bool hasCrossSectionData() const; + + //! Enters the collision data found in the YAML when it is given in the old (still accepted but deprecated) format. + void applyCollisionData(const AnyMap& node); + + //! Checks the validity of the YAML entry. + void validateCollisionData(const AnyMap& node) const; + private: //! The name of the kind of electron collision string m_kind; @@ -222,7 +234,7 @@ class ElectronCollisionPlasmaRate : public ReactionRate string m_product; //! The energy threshold of electron collision - double m_threshold; + double m_threshold = 0; //! electron energy levels [eV] vector m_energyLevels; @@ -240,11 +252,19 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! collision cross sections [m2] after interpolation vector m_crossSectionsInterpolated; - //! collision cross section [m2] interpolated on #m_energyLevels offset by the - //! threshold energy (the first energy level). - //! This is used for the calculation of the super-elastic collision reaction - //! rate coefficient. + //! Super-elastic cross sections [m²] interpolated on the shared EEDF grid + //! after shifting the tabulated energy levels by #m_threshold. Eigen::ArrayXd m_crossSectionsOffset; + + //! The name of the collision field in the new YAML PlasmaPhase implementation + string m_collisionName; + + //! Check whether a YAML node entry offers some cross section data + bool m_hasCrossSectionData = false; + + //! Get the energy threshold for the reaction from the energy levels and cross sections data if threshold is not given in the reaction + void setDefaultThreshold(); + }; } diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 7f17e9de723..b91531c01bc 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -48,12 +48,14 @@ struct TwoTempPlasmaData : public ReactionData * activation energy for electron is included. * * @f[ - * k_f = A T_e^b \exp (-E_{a,g}/RT) \exp (E_{a,e} (T_e - T)/(R T T_e)) + * k_f = A T^{b_g} T_e^b * exp(-E_{a,g}/RT) * exp(E_{a,e}(T_e - T)/(R T T_e)) * @f] * * where @f$ T_e @f$ is the electron temperature, @f$ E_{a,g} @f$ is the activation * energy for gas, and @f$ E_{a,e} @f$ is the activation energy for electron, see * Kossyi, et al. @cite kossyi1992. + * The optional gas temperature exponent b_g defaults to zero, which stricly corresponds to @cite kossyi1992. + * If b_g is non-zero, a generalisation is used. * * @ingroup arrheniusGroup */ @@ -69,7 +71,10 @@ class TwoTempPlasmaRate : public ArrheniusBase * @param b Temperature exponent (non-dimensional) * @param Ea Activation energy in energy units [J/kmol] * @param EE Activation electron energy in energy units [J/kmol] + * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. */ + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg); + TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); @@ -90,12 +95,12 @@ class TwoTempPlasmaRate : public ArrheniusBase */ double evalFromStruct(const TwoTempPlasmaData& shared_data) const { // m_E4_R is the electron activation (in temperature units) - return m_A * std::exp(m_b * shared_data.logTe - - m_Ea_R * shared_data.recipT + - m_E4_R * (shared_data.electronTemp - shared_data.temperature) - * shared_data.recipTe * shared_data.recipT); + return m_A * std::exp(m_bg * shared_data.logT + m_b * shared_data.logTe + - m_Ea_R * shared_data.recipT + m_E4_R * (shared_data.electronTemp - shared_data.temperature) + * shared_data.recipTe * shared_data.recipT); } + //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! @@ -109,6 +114,9 @@ class TwoTempPlasmaRate : public ArrheniusBase double activationElectronEnergy() const { return m_E4_R * GasConstant; } + +protected: + double m_bg = 0.0; //!< Gas temperature exponent }; } diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h new file mode 100644 index 00000000000..2df4def6e8a --- /dev/null +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -0,0 +1,293 @@ +//! @file VibrationalRelaxationRate.h +//! Header for vibrational relaxation reaction rates in plasma kinetics. +//! +//! This file is part of Cantera. See License.txt in the top-level directory or +//! at https://cantera.org/license.txt for license and copyright information. + + +//! It implements the relaxation rate of vibrationally excited species in plasma kinetics. +//! Four models are currently supported by this class: 'constant', 'multi-state-resolved', 'starikovskiy', and 'castela'. + +//! The 'constant' model relaxes the vibrational species with a constant rate coefficient. It could just as well be an +//! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. + +//! The 'multi-state-resolved' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species +//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates are based on the SSH theory +//! detailed in Chapter 7 of the following book: +//! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). +//! Plasma kinetics in atmospheric gases (Vol. 31). Springer Science & Business Media. +//! The simplified version of the SSH theory implemented in this model is based on the harmonic oscillator approximation and can be found +//! in the following paper (equations 18 and 19): +//! Vasco Guerra et al 2019 Plasma Sources Sci. Technol. 28 073001 +//! The k10 rates are taken from //! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted +//! low-temperature NH3 oxidation,” Combustion and Flame, 256, which itself took those rates from both the Capitelli book and +//! the article from Strikovskiy and Aleksandrov (2013), see more below. + +//! The 'castela' model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. +//! It implements the mean vibrational energy relaxation model using a fictive cantera species. +//! The Castela model is based on the following paper: +//! Castela, M., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., & Laux, C. O. (2016). +//! Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of +//! plasma-assisted ignition in turbulent flows. Combustion and flame, 166, 133-147. + +//! The so-called 'starikovskiy' model is an extension of the Castela model to several vibrational species and additional colliders. +//! A lot of vibrational relaxation rates can be found in Table 1 of the following paper: +//! Starikovskiy, A., & Aleksandrov, N. (2013). Plasma-assisted ignition and combustion. +//! Progress in Energy and Combustion Science, 39(1), 61-110. +//! The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material +//! of the following paper: +//! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted low-temperature NH3 oxidation,” +//! Combustion and Flame, 256. +//! More rates for the vibrational relaxation of CH$ can be found in the following paper: +//! Popov, N.A. (2016) “Kinetics of plasma-assisted combustion: Effect of non-equilibrium excitation on the ignition +//! and oxidation of combustible mixtures,” Plasma Sources Science and Technology. Institute of Physics Publishing. + + + +#ifndef CT_VIBRATIONALRELAXATIONRATE_H +#define CT_VIBRATIONALRELAXATIONRATE_H + +#include "Arrhenius.h" + +#include + +namespace Cantera +{ + +//! Shared temperature data for vibrational relaxation rates. +struct DetailedVibData : public ReactionData +{ + //! Update cached temperature-dependent data. + /*! + * @param phase Thermodynamic phase used to retrieve the gas temperature. + * @param kin Kinetics object. Not used here, but required by the + * ReactionData interface. + * + * @return `true` if the temperature has changed and rates need to be + * recomputed; `false` otherwise. + */ + bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; +}; + + +//! Vibrational relaxation reaction rate. +/*! + * This class provides a common implementation for several vibrational + * relaxation models: + * + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + * + * Internally, all models are mapped to the following generic expression: + * + * @f[ + * k_f = + * scaling \, A \, + * \exp \left( + * b \ln T + * + B + * + C T^{-1/3} + * + D T^{-m} + * + E T^{-z} + * \right) + * @f] + * + * where `T` is the gas temperature in K. + * + * The YAML reaction type is: + * + * @code{.yaml} + * type: vibrational-relaxation + * @endcode + * + * The selected physical model is specified separately using: + * + * @code{.yaml} + * vibration_model: constant + * vibration_model: multi-state-resolved + * vibration_model: starikovskiy + * vibration_model: castela + * @endcode + * + * Unit conventions: + * + * - `A` uses standard Cantera rate coefficient units. Its units depend on the + * reaction order and are converted by `ArrheniusBase`. + * - `b`, `B`, `m`, `z`, and `scaling` are dimensionless. + * - `C` is interpreted as K^(1/3), assuming `T` is in K. + * - `D` is interpreted as K^m, assuming `T` is in K. + * - `E` is interpreted as K^z, assuming `T` is in K. + * + * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as + * raw floating-point values. They are not converted by Cantera's unit system. + * + * @ingroup arrheniusGroup + */ +class VibrationalRelaxationRate : public ArrheniusBase +{ +public: + //! Default constructor. + VibrationalRelaxationRate(); + + //! Constructor using the internal generic representation. + /*! + * @param A Pre-exponential factor. + * @param B Dimensionless constant in the exponential. + * @param C Coefficient multiplying T^(-1/3). + * @param D Coefficient multiplying T^(-m). + * @param b Dimensionless temperature exponent. + * @param scaling Dimensionless scaling factor. + * @param m Temperature exponent used by the D term. + * @param E Coefficient multiplying T^(-z). + * @param z Temperature exponent used by the E term. + */ + VibrationalRelaxationRate(double A, double B, double C, double D, + double b, double scaling = 1.0, + double m = 2.0 / 3.0, + double E = 0.0, double z = 1.0); + + //! Constructor based on AnyMap content. + explicit VibrationalRelaxationRate(const AnyMap& node, + const UnitStack& rate_units = {}); + + //! Set rate parameters from an AnyMap. + void setParameters(const AnyMap& node, + const UnitStack& rate_units) override; + + //! Get rate parameters for YAML serialization. + void getParameters(AnyMap& node) const override; + + //! Create a rate evaluator for this reaction rate type. + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! String identifying this reaction rate type. + const string type() const override { + return "vibrational-relaxation"; + } + + //! Set context of reaction rate evaluation. + /*! + * Vibrational relaxation rates are intended for irreversible + * non-equilibrium plasma reactions. Reversible reactions are rejected + * because the reverse rate cannot be obtained from conventional + * thermochemistry for these models. + */ + void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate the forward rate coefficient. + double evalFromStruct(const DetailedVibData& shared_data) const { + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_scaling * m_A * std::exp( + m_b * shared_data.logT + + m_B + + m_C * invT13 + + m_D * std::pow(invT, m_m) + + m_E * std::pow(invT, m_z) + ); + } + + //! Evaluate the scaled temperature derivative. + /*! + * This returns: + * + * @f[ + * \frac{1}{k_f} \frac{d k_f}{dT} + * = + * \frac{d \ln k_f}{dT} + * @f] + * + * For the internal generic expression, this is: + * + * @f[ + * \frac{b}{T} + * - \frac{C}{3} T^{-4/3} + * - m D T^{-m-1} + * - z E T^{-z-1} + * @f] + */ + double ddTScaledFromStruct(const DetailedVibData& shared_data) const; + +private: + //! Dimensionless constant in the exponential. + double m_B = 0.0; + + //! Coefficient multiplying T^(-1/3). + double m_C = 0.0; + + //! Coefficient multiplying T^(-m). + double m_D = 0.0; + + //! Dimensionless scaling factor. + double m_scaling = 1.0; + + //! Temperature exponent used by the D term. + double m_m = 2.0 / 3.0; + + //! Coefficient multiplying T^(-z). + double m_E = 0.0; + + //! Temperature exponent used by the E term. + double m_z = 1.0; + + //! Castela coefficient a. + double m_castela_a = 0.0; + + //! Castela coefficient b. + double m_castela_b = 0.0; + + //! Castela reference pressure. + double m_referencePressure = OneAtm; + + //! Selected vibrational relaxation model. + /*! + * Accepted values: + * + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + */ + string m_vibration_model = "multi-state-resolved"; + + //! YAML key names. + string m_B_str = "B"; + string m_C_str = "C"; + string m_D_str = "D"; + string m_scaling_str = "scaling"; + string m_m_str = "m"; + string m_E_str = "E"; + string m_z_str = "z"; + string m_reference_pressure_str = "reference-pressure"; + + //! Configure the ArrheniusBase part from an already-converted internal A value. + /*! + * This is needed for models such as Castela, where the user-facing YAML does + * not contain a standard Arrhenius A coefficient. + */ + void configureBaseFromInternalA(const AnyMap& node, + const UnitStack& rate_units, + double A, double b); + + //! Configure the ArrheniusBase part from a YAML A value and an explicit b. + /*! + * This is needed for models such as constant and Starikovskiy, where the YAML + * does not contain the standard Arrhenius pair A / b. + */ + void configureBaseFromYamlA(const AnyMap& node, + const UnitStack& rate_units, + const AnyValue& A, + double b); + + }; + +} + +#endif \ No newline at end of file diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index fddc47d5448..b22fa500584 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -61,28 +61,119 @@ class EEDFTwoTermApproximation virtual ~EEDFTwoTermApproximation() = default; - // compute the EEDF given an electric field - // CQM The solver will take the species to consider and the set of cross-sections - // from the PlasmaPhase object. - // It will write the EEDF and its grid into the PlasmaPhase object. - // Successful returns are indicated by a return value of 0. + //! Compute the electron energy distribution function for the current plasma state. + /*! + * The solver uses the electric field, species mole fractions, and electron + * collision cross sections provided by the associated PlasmaPhase object. + * On success, the internal EEDF and energy grid are updated. + * + * @return 0 if the EEDF calculation succeeds. + */ int calculateDistributionFunction(); + //! Sets a linear energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. void setLinearGrid(double& kTe_max, size_t& ncell); + //! Sets a quadratic energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. + void setQuadraticGrid(double& kTe_max, size_t& ncell); + + //! Sets a geometric energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. + //! the geometric ratio defaults to 1.01 if not specified. + void setGeometricGrid(double& kTe_max, size_t& ncell, double ratio = 1.01); + + //! Sets a custom energy grid for the EEDF solver, defined by the user-provided vector of energy levels. + void setCustomGrid(span levels); + + //! Build or rebuild the grid-dependent cache used for scattering matrices. void setGridCache(); + //! Set the type of generated electron energy grid. + /*! + * Supported values are "Linear", "Quadratic", and "Geometric". + * + * @param gridType Type of grid spacing to use when generating or adapting + * the electron energy grid. + */ + void setGridType(const string& gridType); + + //! Set the initial grid parameters used by generated EEDF grids. + /*! + * These values are used when creating Linear, Quadratic, or Geometric grids, + * and are also reused during grid adaptation. + * + * @param initialMaxEnergy Maximum electron energy of the initial grid [eV]. + * @param nGridCells Number of grid cells. The number of grid edges is + * nGridCells + 1. + */ + void setInitialGridParameters(double initialMaxEnergy, size_t nGridCells); + + //! This funcion enables or disables the grid adaptation based on the EEDF decay at its tail. + void enableGridAdaptation(bool enabled); + + //! Set parameters controlling automatic adaptation of the EEDF energy grid. + /*! + * Grid adaptation adjusts the maximum grid energy based on the number of + * decades by which the EEDF decays between the low- and high-energy ends of + * the grid. The number of grid cells is kept fixed. + * + * @param enabled True to enable grid adaptation. + * @param minDecayDecades Minimum acceptable EEDF tail decay in decades. If + * the decay is smaller, the maximum grid energy is + * increased. + * @param maxDecayDecades Maximum acceptable EEDF tail decay in decades. If + * the decay is larger, the maximum grid energy is + * decreased. + * @param updateFactor Relative factor used to increase or decrease the + * maximum grid energy during adaptation. + * @param maxIterations Maximum number of grid adaptation iterations per + * EEDF solve. + */ + void setGridAdaptationParameters(bool enabled, + double minDecayDecades, + double maxDecayDecades, + double updateFactor, + size_t maxIterations); + + //! Return the electron energy grid edges [eV]. + /*! + * The returned vector contains m_points + 1 values corresponding to cell + * boundaries. + */ vector getGridEdge() const { return m_gridEdge; } + //! Return the EEDF values interpolated at the electron energy grid edges. + /*! + * These values are copied back to PlasmaPhase after a successful + * Boltzmann-two-term EEDF solve. + */ vector getEEDFEdge() const { return m_f0_edge; } - double getElectronMobility() const { - return m_electronMobility; - } + //! Return the latest value of the computed electron mobility computed from the EEDF + double getElectronMobility() const; + + //! Controls the threshold below which the EEDF solver does not solve for the EEDF but imposes a Maxwellian distribution. + void setReducedFieldThresholdBeforeMaxwellianTd(double threshold); + + //! An extension of the linearInterp function that returns specified values when the input is + //! out of bounds instead of returning one of the extremities of the list. + double linearInterpBounded(double x, + span xpts, + span fpts, + double below_value, + double above_value); + + //! A special use case of linearInterpBounded for cross sections + //! It is not necessary as the code could just use the previous function + //! using the values 0, but the author of this code finds it clearer to + //! the reader to have an explicit name for this use case. + double linearInterpCrossSectionZeroOutside(double x, + span xpts, + span fpts); + protected: @@ -192,7 +283,7 @@ class EEDFTwoTermApproximation double electronDiffusivity(const Eigen::VectorXd& f0); //! Mobility - double electronMobility(const Eigen::VectorXd& f0); + double computeElectronMobility(const Eigen::VectorXd& f0); //! Initialize species indices associated with cross-section data void initSpeciesIndexCrossSections(); @@ -209,13 +300,16 @@ class EEDFTwoTermApproximation //! Compute the total (elastic + inelastic) cross section void calculateTotalCrossSection(); + //! + void checkSpeciesNoCrossSection(); + //! Compute the L1 norm of a function f defined over a given energy grid. //! //! @param f Vector representing the function values (EEDF) //! @param grid Vector representing the energy grid corresponding to f double norm(const Eigen::VectorXd& f, const Eigen::VectorXd& grid); - //! Electron mobility [m²/V·s] + //! Electron mobility computed from the most recent valid EEDF [m²/V/s]. double m_electronMobility; //! Grid of electron energy (cell center) [eV] @@ -236,7 +330,7 @@ class EEDFTwoTermApproximation //! The energy boundaries of the overlap of cell i and j vector>> m_eps; - //! normalized electron energy distribution function + //! Normalized electron energy distribution function Eigen::VectorXd m_f0; //! EEDF at grid edges (cell boundaries) @@ -249,13 +343,13 @@ class EEDFTwoTermApproximation //! energy grid vector m_totalCrossSectionEdge; - //! vector of total elastic cross section weighted with mass ratio + //! Vector of total elastic cross section weighted with mass ratio vector m_sigmaElastic; - //! list of target species indices in global Cantera numbering (1 index per cs) + //! List of target species indices in global Cantera numbering (1 index per cs) vector m_kTargets; - //! list of target species indices in local X EEDF numbering (1 index per cs) + //! List of target species indices in local X EEDF numbering (1 index per cs) vector m_klocTargets; //! Indices of species which has no cross-section data @@ -270,10 +364,14 @@ class EEDFTwoTermApproximation //! Previous mole fraction of targets used to compute eedf vector m_X_targets_prev; - //! in factor. This is used for calculating the Q matrix of - //! scattering-in processes. + //! Multiplicity factor for scattering-in terms in matrix_Q(). + /*! + * Values are 2 for ionization, 0 for attachment, and 1 for other collision + * types. + */ vector m_inFactor; + //! Conversion factor sqrt(2 e / m_e) used in electron energy-space transport terms. double m_gamma; //! flag of having an EEDF @@ -281,6 +379,43 @@ class EEDFTwoTermApproximation //! First call to calculateDistributionFunction bool m_first_call; + + //! The grid type used by the EEDF solver. Supported values are "Linear", "Quadratic", and "Geometric". + string m_gridType = "Linear"; + + //! Initial maximum energy of the grid [eV] used when generating EEDF grids. Can be changed by user input of by grid adaptation. + double m_kTeMax = 60.0; + + //! Initial / default number of grid cells. Can be changed by user input. + size_t m_initialGridCells = 301; + + //! Flag indicating whether grid adaptation based on EEDF tail decay is enabled. + bool m_adaptGrid = false; + + // Parameters controlling grid adaptation based on EEDF tail decay. They will be used if no parameters are provided by the user in the YAML file. + + //! Minimum acceptable EEDF tail decay in decades; lower values trigger grid expansion. + double m_minEedfDecay = 10.0; + + //! Maximum acceptable EEDF tail decay in decades; higher values trigger grid contraction. + double m_maxEedfDecay = 12.0; + + //! Relative factor used to increase or decrease the maximum grid energy. + double m_gridUpdateFactor = 0.25; + + //! Maximum number of grid adaptation iterations per EEDF solve. + size_t m_maxGridAdaptIterations = 100; + + //! Sets a new grid for the EEDF solver once the maximum energy is updated. This function is only called in the case of grid adaptation. + void updateGrid(double maxEnergy); + + //! Reduced electric field threshold below which the EEDF is forced to a Maxwellian. + /*! + * Stored internally in SI units [V m²]. User-facing input is provided in + * Townsend by setReducedFieldThresholdBeforeMaxwellianTd(). + */ + double EN_min = 1e-21; + }; // end of class EEDFTwoTermApproximation } // end of namespace Cantera diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 61bb85c9338..2017356e7b3 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -22,7 +22,7 @@ class ElectronCollisionPlasmaRate; //! Base class for handling plasma properties, specifically focusing on the //! electron energy distribution. /*! - * This class provides functionality to manage the the electron energy distribution + * This class provides functionality to manage the electron energy distribution * using two primary methods for defining the electron distribution and electron * temperature. * @@ -113,7 +113,10 @@ class PlasmaPhase: public IdealGasPhase //! Get electron energy levels. //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints void getElectronEnergyLevels(span levels) const { - Eigen::Map(levels.data(), levels.size()) = m_electronEnergyLevels; + checkArraySize("PlasmaPhase::getElectronEnergyLevels", + levels.size(), m_nPoints); + Eigen::Map(levels.data(), levels.size()) = + m_electronEnergyLevels; } //! Set discretized electron energy distribution. @@ -128,7 +131,10 @@ class PlasmaPhase: public IdealGasPhase //! @param distrb The vector of electron energy distribution. //! Length: #m_nPoints. void getElectronEnergyDistribution(span distrb) const { - Eigen::Map(distrb.data(), distrb.size()) = m_electronEnergyDist; + checkArraySize("PlasmaPhase::getElectronEnergyDistribution", + distrb.size(), m_nPoints); + Eigen::Map(distrb.data(), distrb.size()) = + m_electronEnergyDist; } //! Set the shape factor of isotropic electron energy distribution. @@ -186,6 +192,7 @@ class PlasmaPhase: public IdealGasPhase return m_do_normalizeElectronEnergyDist; } + //! adds a species to the phase. Override from IdealGasPhase to take into account electrons. bool addSpecies(shared_ptr spec) override; //! Electron Temperature (K) @@ -337,7 +344,7 @@ class PlasmaPhase: public IdealGasPhase return m_levelNum; } - //! Get the indicies for inelastic electron collisions + //! Get the indices for inelastic electron collisions //! @since New in %Cantera 3.2. const vector& kInelastic() const { return m_kInelastic; @@ -368,10 +375,15 @@ class PlasmaPhase: public IdealGasPhase //! Set the absolute electric field strength [V/m] void setElectricField(double E) { + if (!std::isfinite(E) || E < 0.0) { + throw CanteraError("PlasmaPhase::setElectricField", + "Electric field must be finite and non-negative."); + } m_electricField = E; } - //! Calculate the degree of ionization + // @todo Add the method to compute the degree of ionization of the plasma. + //!Calculate the degree of ionization //double ionDegree() const { // double ne = concentration(m_electronSpeciesIndex); // [kmol/m³] // double n_total = molarDensity(); // [kmol/m³] @@ -385,7 +397,18 @@ class PlasmaPhase: public IdealGasPhase //! Set reduced electric field given in [V·m²] void setReducedElectricField(double EN) { - m_electricField = EN * molarDensity() * Avogadro; // [V/m] + if (!std::isfinite(EN) || EN < 0.0) { + throw CanteraError("PlasmaPhase::setReducedElectricField", + "Reduced electric field must be finite and non-negative."); + } + + const double nDensity = molarDensity() * Avogadro; + if (!std::isfinite(nDensity) || nDensity <= 0.0) { + throw CanteraError("PlasmaPhase::setReducedElectricField", + "Cannot set reduced electric field with non-positive number density."); + } + + m_electricField = EN * nDensity; } virtual void setSolution(std::weak_ptr soln) override; @@ -428,8 +451,36 @@ class PlasmaPhase: public IdealGasPhase void endEquilibrate() override; + //! Calculates the intrinsic heating power (W/m³) of the plasma. double intrinsicHeating() override; + //! Calculates the power losses (W/m³) of the plasma electrons through inelastic collisions. + double inelasticPower(); + + //! Allows to store info on potential vibrational reservoir / fictive species that could be declare by the user + //! These species take the form A(v) where A is a base species and v represents the general vibration. + struct VibrationalReservoirSpecies { + size_t reservoirIndex = npos; + size_t baseSpeciesIndex = npos; + }; + + //! list of detected vibrational reservoir species in the phase. + std::vector m_vibrationalReservoirSpecies; + + //! Flag to indicate whether the vibrational reservoir species list needs to be updated + bool m_vibrationalReservoirSpeciesNeedUpdate = true; + + // Thresholds for identifying vibrational reservoir species based on their mole fractions + double m_vibrationalMoleFractionThreshold = 1e-2; + double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + + //! Update the list of vibrational reservoir / fictive species + void updateVibrationalReservoirSpecies(); + + //! Checks that declared vibrational reservoir / fictive species have mole fractions + //! low enough not to impact the phase chemistry. + void checkVibrationalReservoirMoleFractions(); + protected: void updateThermo() const override; @@ -438,7 +489,7 @@ class PlasmaPhase: public IdealGasPhase void electronEnergyDistributionChanged(); //! When electron energy level changed, plasma properties such as - //! electron-collision reaction rates need to be re-evaluate. + //! electron-collision reaction rates need to be re-evaluated. //! In addition, the cross-sections need to be interpolated at //! the new level. void electronEnergyLevelChanged(); @@ -452,7 +503,7 @@ class PlasmaPhase: public IdealGasPhase //! Check the electron energy distribution /*! - * This method check the electron energy distribution for the criteria + * This method checks the electron energy distribution for the criteria * below. * * 1. The values of electron energy distribution cannot be negative. @@ -527,6 +578,7 @@ class PlasmaPhase: public IdealGasPhase //! where i is the specific process, j is the index of vector. Unit: [eV] vector> m_energyLevels; + // @todo Add a variable to track the ionization degree of the plasma. //! ionization degree for the electron-electron collisions (tmp is the previous one) //double m_ionDegree = 0.0; @@ -560,6 +612,11 @@ class PlasmaPhase: public IdealGasPhase */ void updateElasticElectronEnergyLossCoefficients(); + //! Add an electron collision with respect to the old formats of electron collision. Function made for compatibility. + void addStandaloneElectronCollision(const AnyMap& item); + + string inferElectronCollisionKind(const shared_ptr& collision) const; + private: //! Solver used to calculate the EEDF based on electron collision rates @@ -600,6 +657,13 @@ class PlasmaPhase: public IdealGasPhase //! Work array mutable std::vector m_work; + + //! The array containing the names of the electron collisions + map m_electronCollisionDefinitions; + + //! The array containing the references of each electron collision. + set m_referencedElectronCollisions; + }; } diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index c80bed26316..53eb23113e2 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -213,6 +213,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double electricField() void updateElectronEnergyDistribution() except +translate_exception double elasticPowerLoss() except +translate_exception + double electronMobility() cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 021f80ca359..b4529abf97f 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1991,6 +1991,17 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) return self.plasma.elasticPowerLoss() + property electron_mobility: + """ + Electron mobility (m^2/(V.s)) + + .. versionadded:: 4.0 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electronMobility() + cdef class InterfacePhase(ThermoPhase): """ A class representing a surface, edge phase """ def __cinit__(self, *args, **kwargs): diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index ce4a1ecbcd6..a21f4fd5e0f 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -31,12 +31,12 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic return false; } - // Update distribution + // Update the cached eedf from the plasma phase. m_dist_number = pp.distributionNumber(); distribution.resize(pp.nElectronEnergyLevels()); pp.getElectronEnergyDistribution(distribution); - // Update energy levels + // Update the cached energy grid only when the phase grid revision changes. if (pp.levelNumber() != levelNumber || energyLevels.empty()) { levelNumber = pp.levelNumber(); energyLevels.resize(pp.nElectronEnergyLevels()); @@ -45,35 +45,56 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic return true; } -void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) +void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) { + ReactionRate::setParameters(node, rate_units); - if (!node.hasKey("energy-levels") && !node.hasKey("cross-sections")) { - return; - } - if (node.hasKey("kind")) { - m_kind = node["kind"].asString(); - } - if (node.hasKey("target")) { - m_target = node["target"].asString(); + if (node.hasKey("collision")) { + m_collisionName = node["collision"].asString(); } - if (node.hasKey("product")) { - m_product = node["product"].asString(); + + bool hasInlineCrossSectionData = + node.hasKey("energy-levels") && node.hasKey("cross-sections"); + + bool isNamedElectronCollision = node.hasKey("name"); + bool isCollisionReference = node.hasKey("collision"); + + if (hasInlineCrossSectionData) { + if (!isNamedElectronCollision && !isCollisionReference) { + writelog("CAREFUL! Inline electron-collision cross-section data without a " + "'name' entry are deprecated. Please move these data to a named " + "'electron-collisions' entry and reference it using 'collision'.\n"); + } + + applyCollisionData(node); } - m_energyLevels = node["energy-levels"].asVector(); - m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); - m_threshold = node.getDouble("threshold", 0.0); } -void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { +void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const +{ node["type"] = type(); + node["energy-levels"] = m_energyLevels; node["cross-sections"] = m_crossSections; + + if (m_threshold != 0.0) { + node["threshold"] = m_threshold; + } + if (!m_kind.empty()) { node["kind"] = m_kind; } + + if (!m_target.empty()) { + node["target"] = m_target; + } + + if (!m_product.empty()) { + node["product"] = m_product; + } } void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( @@ -142,7 +163,7 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( m_crossSectionsOffset.resize(shared_data.energyLevels.size()); for (size_t i = 1; i < m_energyLevels.size(); i++) { // The energy levels are offset by the first energy level (threshold) - superElasticEnergyLevels.push_back(m_energyLevels[i] - m_energyLevels[0]); + superElasticEnergyLevels.push_back(m_energyLevels[i] - m_threshold); } for (size_t i = 0; i < shared_data.energyLevels.size(); i++) { // The interpolated super-elastic cross section is evaluated @@ -159,19 +180,20 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( shared_data.energyLevels.data(), shared_data.energyLevels.size() ); - // Map energyLevels in Eigen::ArrayXd + // Map the electron energy distribution to Eigen::ArrayXd. auto distribution = Eigen::Map( shared_data.distribution.data(), shared_data.distribution.size() ); // unit in kmol/m3/s kr = pow(2.0 * ElectronCharge / ElectronMass, 0.5) * Avogadro * - simpson((eps + m_energyLevels[0]).cwiseProduct( + simpson((eps + m_threshold).cwiseProduct( distribution.cwiseProduct(m_crossSectionsOffset)), eps); } void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) { + const ThermoPhase& thermo = kin.thermo(); // get electron species name string electronName; @@ -215,16 +237,7 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } - if (m_threshold == 0.0 && - (m_kind == "excitation" || m_kind == "ionization" || m_kind == "attachment")) - { - for (size_t i = 0; i < m_energyLevels.size(); i++) { - if (m_energyLevels[i] > 0.0) { // Look for first non-zero cross-section - m_threshold = m_energyLevels[i]; - break; - } - } - } + setDefaultThreshold(); if (!rxn.reversible) { return; // end checking of forward reaction @@ -245,4 +258,118 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } +void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) +{ + if (node.hasKey("kind")) { + string collisionKind = node["kind"].asString(); + + if (!m_kind.empty() && m_kind != collisionKind) { + string collisionName = node.hasKey("name") ? node["name"].asString() : m_collisionName; + + bool allowCollapsedInelastic = + (m_kind == "effective" || m_kind == "elastic") && + collisionKind == "excitation"; + + if (!allowCollapsedInelastic) { + throw InputFileError("applyCollisionData", node, + "Electron collision '{}' has kind '{}', but the reaction was inferred as '{}'.", + collisionName, collisionKind, m_kind); + } + + warn_user("ElectronCollisionPlasmaRate::applyCollisionData", + "Electron collision '{}' has kind '{}', but the reaction was inferred " + "as '{}'. Treating this as an intentionally collapsed inelastic " + "electron-collision channel. No species source term will be generated " + "for the unresolved product.", + collisionName, collisionKind, m_kind); + } + + // Important: keep the collision classified using the data-file kind. + // For collapsed channels such as N2(rot), this ensures that the + // Boltzmann solver treats the channel as inelastic. + m_kind = collisionKind; + } + + if (node.hasKey("target")) { + m_target = node["target"].asString(); + } + + if (node.hasKey("product")) { + m_product = node["product"].asString(); + } + + if (!node.hasKey("energy-levels")) { + throw InputFileError("applyCollisionData", node, "Missing 'energy-levels'"); + } + + if (!node.hasKey("cross-sections")) { + throw InputFileError("applyCollisionData", node, "Missing 'cross-sections'"); + } + + m_energyLevels = node["energy-levels"].asVector(); + m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); + m_threshold = node.getDouble("threshold", 0.0); + + setDefaultThreshold(); + + validateCollisionData(node); + m_hasCrossSectionData = true; +} + +void ElectronCollisionPlasmaRate::validateCollisionData(const AnyMap& node) const +{ + if (m_energyLevels.size() < 2) { + throw InputFileError("validateCollisionData" , node, "Need at least two energy levels."); + } + + if (m_energyLevels.size() != m_crossSections.size()) { + throw InputFileError("validateCollisionData" , node, "energy-levels and cross-sections size mismatch."); + } + + for (size_t i = 0; i < m_energyLevels.size(); i++) { + if (!std::isfinite(m_energyLevels[i]) || m_energyLevels[i] < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative energy level value"); + } + if (!std::isfinite(m_crossSections[i]) || m_crossSections[i] < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative cross-section value"); + } + if (i > 0 && m_energyLevels[i] <= m_energyLevels[i - 1]) { + throw InputFileError("validateCollisionData" , node, "energy-levels must be strictly increasing."); + } + } + + if (!std::isfinite(m_threshold) || m_threshold < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative threshold value"); + } +} + +void ElectronCollisionPlasmaRate::setDefaultThreshold() +{ + if (m_threshold != 0.0 || m_energyLevels.empty()) { + return; + } + + if (m_kind != "excitation" && m_kind != "ionization" && m_kind != "attachment") { + return; + } + + for (double level : m_energyLevels) { + if (level > 0.0) { + m_threshold = level; + break; + } + } +} + +const string& ElectronCollisionPlasmaRate::collisionName() const +{ + return m_collisionName; +} + +bool ElectronCollisionPlasmaRate::hasCrossSectionData() const +{ + return m_hasCrossSectionData; +} + + } diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 1019fe4f4ed..cb079dc0ee1 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,6 +18,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" namespace Cantera { @@ -40,6 +41,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // VibrationalRelaxationRate evaluator + reg("vibrational-relaxation", [](const AnyMap& node, const UnitStack& rate_units) { + return new VibrationalRelaxationRate(node, rate_units); + }); + // ElectronCollisionPlasmaRate evaluator reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) { return new ElectronCollisionPlasmaRate(node, rate_units); diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index ce3f4694969..0b6563e98e3 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -57,19 +57,34 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; m_E4_R = EE / GasConstant; + m_bg = 0.0; +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) + : ArrheniusBase(A, b, Ea) +{ + TwoTempPlasmaRate(A, b, Ea, EE); + m_bg = bg; } TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units) : TwoTempPlasmaRate() { setParameters(node, rate_units); + + if (node.hasKey("b-gas")) { + m_bg = node["b-gas"].asDouble(); + } else if (node.hasKey("b_gas")) { + m_bg = node["b_gas"].asDouble(); + } } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); - return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; + return m_bg * shared_data.recipT + + (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp new file mode 100644 index 00000000000..79aabf911ee --- /dev/null +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -0,0 +1,835 @@ +//! @file VibrationalRelaxationRate.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/VibrationalRelaxationRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/global.h" + +#include +#include +#include +#include +#include +#include + +namespace Cantera +{ + +namespace +{ + +void requireNoKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (node.hasKey(key)) { + throw InputFileError(where, node, + "Key '{}' is not allowed for vibration_model '{}'.", key, model); + } +} + + +void requireKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (!node.hasKey(key)) { + throw InputFileError(where, node, + "Missing required key '{}' for vibration_model '{}'.", key, model); + } +} + + +bool isVibrationalSpecies(const string& name) +{ + // Supported names: + // + // N2(v) + // O2(v0) + // O2(v1) + // O2(v12) + // + // The old form O2(v=1) is intentionally not required. + const auto pos = name.find("(v"); + return pos != string::npos && !name.empty() && name.back() == ')'; +} + + +string groundStateName(const string& name) +{ + const auto pos = name.find("(v"); + if (pos == string::npos) { + return name; + } + return name.substr(0, pos); +} + + +string vibrationalFamilyName(const string& name) +{ + // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v). + return groundStateName(name) + "(v)"; +} + + +double compositionSum(const Composition& comp) +{ + double sum = 0.0; + for (const auto& item : comp) { + sum += item.second; + } + return sum; +} + + +Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) +{ + Composition out; + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (isVibrationalSpecies(name)) { + out[groundStateName(name)] += value; + } else { + out[name] += value; + } + } + return out; +} + + +bool sameComposition(const Composition& a, const Composition& b, + double tol = 1e-12) +{ + std::set names; + + for (const auto& item : a) { + names.insert(item.first); + } + for (const auto& item : b) { + names.insert(item.first); + } + + for (const auto& name : names) { + double av = 0.0; + double bv = 0.0; + + const auto ait = a.find(name); + if (ait != a.end()) { + av = ait->second; + } + + const auto bit = b.find(name); + if (bit != b.end()) { + bv = bit->second; + } + + if (std::abs(av - bv) > tol) { + return false; + } + } + + return true; +} + + +std::vector vibrationalSpeciesInComposition(const Composition& comp) +{ + std::vector out; + + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (value != 0.0 && isVibrationalSpecies(name)) { + out.push_back(name); + } + } + + return out; +} + + +// Registry used to check that a given vibrational family uses exactly one +// relaxation model inside one Kinetics object. +// +// Allowed: +// +// N2(v) -> constant +// O2(v) -> multi-state-resolved +// NH3(v) -> starikovskiy +// +// Forbidden: +// +// N2(v) + O -> castela +// N2(v) + N2 -> starikovskiy +// +std::map> s_modelByFamily; +std::set s_warnedMixedModels; + + +void registerVibrationalModelConsistency(const Kinetics& kin, + const string& family, + const string& model, + const AnyMap& input) +{ + auto& modelByFamily = s_modelByFamily[&kin]; + + const auto existing = modelByFamily.find(family); + if (existing != modelByFamily.end() && existing->second != model) { + throw InputFileError("VibrationalRelaxationRate::setContext", input, + "Inconsistent vibration_model for vibrational family '{}'. " + "This family was already registered with model '{}', but the " + "current reaction uses model '{}'. A given vibrational family " + "must use exactly one relaxation model.", + family, existing->second, model); + } + + modelByFamily[family] = model; + + std::set models; + for (const auto& item : modelByFamily) { + models.insert(item.second); + } + + if (models.size() > 1 && !s_warnedMixedModels.count(&kin)) { + std::ostringstream msg; + msg << "Multiple vibrational relaxation models were detected in the " + << "same kinetics object. This is allowed only if each " + << "vibrational family is internally consistent:\n"; + + for (const auto& item : modelByFamily) { + msg << " - " << item.first << ": " << item.second << "\n"; + } + + warn_user("VibrationalRelaxationRate::setContext", msg.str()); + s_warnedMixedModels.insert(&kin); + } +} + + +string inferRelaxingFamily(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A vibrational-relaxation reaction must contain at least one " + "vibrational reactant, for example O2(v1), N2(v), or NH3(v1)."); + } + + return vibrationalFamilyName(vibReactants.front()); +} + + +void validateSimpleRelaxationToGroundState(const Reaction& rxn, + const string& model) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.size() != 1) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects exactly one vibrational reactant.", + model); + } + + if (!vibProducts.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' describes relaxation to the ground state " + "and therefore does not allow vibrational products.", model); + } + + const Composition relaxedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + if (!sameComposition(relaxedReactants, rxn.products)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid vibrational relaxation stoichiometry for model '{}'. " + "Expected a reaction equivalent to X(v) + M => X + M, where the " + "collider M is unchanged.", model); + } + + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects a bimolecular relaxation reaction " + "of the form X(v) + M => X + M.", model); + } +} + + +void validateCastelaReaction(const Reaction& rxn) +{ + validateSimpleRelaxationToGroundState(rxn, "castela"); + + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const string family = vibrationalFamilyName(vibReactants.front()); + + if (family != "N2(v)") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model is only valid for N2 vibrational " + "relaxation. Found vibrational family '{}'.", family); + } + + // Castela parameters used here are only intended for the following + // colliders: N2, O2, and O. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + string collider = ""; + for (const auto& item : collapsedReactants) { + const string& name = item.first; + const double value = item.second; + + if (name == "N2") { + if (std::abs(value - 2.0) < 1e-12) { + collider = "N2"; + } + } else if (std::abs(value - 1.0) < 1e-12) { + collider = name; + } + } + + if (collider != "N2" && collider != "O2" && collider != "O") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model only supports colliders 'N2', " + "'O2', and 'O'. Found collider '{}'.", collider); + } +} + + +void validateDetailedRelaxationReaction(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A detailed vibrational relaxation reaction must contain at least " + "one vibrational reactant."); + } + + // All vibrational species in a detailed relaxation reaction are required + // to belong to the same vibrational family. + const string family = vibrationalFamilyName(vibReactants.front()); + + for (const auto& sp : vibReactants) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational reactants must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + for (const auto& sp : vibProducts) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational products must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + // The reaction must conserve the ground-state composition once all + // vibrational labels are removed. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + const Composition collapsedProducts = + replaceVibrationalSpeciesByGroundState(rxn.products); + + if (!sameComposition(collapsedReactants, collapsedProducts)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: replacing all " + "vibrational species by their ground-state species does not " + "conserve stoichiometry."); + } + + // The current detailed VV/VT formulation is bimolecular. + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: expected a " + "bimolecular reaction with two reactant molecules and two product " + "molecules."); + } +} + +} // namespace + + +bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + const double T = phase.temperature(); + + if (T == temperature) { + return false; + } + + ReactionData::update(T); + return true; +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate() +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + double A, double B, double C, double D, + double b, double scaling, double m, double E, double z) + : ArrheniusBase(A, b, 0.0) + , m_B(B) + , m_C(C) + , m_D(D) + , m_scaling(scaling) + , m_m(m) + , m_E(E) + , m_z(z) +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + const AnyMap& node, const UnitStack& rate_units) + : VibrationalRelaxationRate() +{ + setParameters(node, rate_units); +} + + +void VibrationalRelaxationRate::configureBaseFromInternalA( + const AnyMap& node, const UnitStack& rate_units, double A, double b) +{ + // Store the original input and configure reaction-rate units. + // + // We intentionally do not call ArrheniusBase::setParameters here because + // some vibration models do not expose a standard YAML rate-constant with + // both A and b. Castela is the main example. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + m_A = A; + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::configureBaseFromYamlA( + const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b) +{ + // Store the original input and configure reaction-rate units first, so + // conversionUnits() is available for A. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + // Convert the user-facing YAML A value with Cantera's standard + // rate-coefficient unit conversion. + m_A = node.units().convertRateCoeff(A, conversionUnits()); + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) +{ + // The reaction rate type is always: + // + // type: vibrational-relaxation + // + // The physical model is selected separately by: + // + // vibration_model: constant + // vibration_model: multi-state-resolved + // vibration_model: starikovskiy + // vibration_model: castela + // + // For backward compatibility, the default model is multi-state-resolved. + m_vibration_model = "multi-state-resolved"; + + if (node.hasKey("vibration_model")) { + m_vibration_model = node["vibration_model"].asString(); + } else if (node.hasKey("model")) { + m_vibration_model = node["model"].asString(); + } + + if (!node.hasKey("rate-constant")) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "A vibrational-relaxation reaction requires a 'rate-constant' " + "mapping."); + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The 'rate-constant' field must be a mapping."); + } + + const auto& rateMap = rate.as(); + + if (m_vibration_model == "constant") { + // Constant model: + // + // k(T) = A + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); + + m_B = 0.0; + m_C = 0.0; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else if (m_vibration_model == "multi-state-resolved") { + // Detailed VV/VT model: + // + // k(T) = scaling * A * exp( + // b * log(T) + // + B + // + C * T^(-1/3) + // + D * T^(-2/3) + // ) + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + ArrheniusBase::setParameters(node, rate_units); + + m_B = rateMap.hasKey(m_B_str) ? rateMap[m_B_str].asDouble() : 0.0; + m_C = rateMap.hasKey(m_C_str) ? rateMap[m_C_str].asDouble() : 0.0; + m_D = rateMap.hasKey(m_D_str) ? rateMap[m_D_str].asDouble() : 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = rateMap.hasKey(m_scaling_str) + ? rateMap[m_scaling_str].asDouble() + : 1.0; + } + else if (m_vibration_model == "starikovskiy") { + // User-facing formula: + // + // k(T) = A * T^n * exp( + // K + // - B * T^(-1/3) + // + C * T^(-m) + // + D * T^(-z) + // ) + // + // Internal mapping: + // + // m_b = n + // m_B = K + // m_C = -B + // m_D = C + // m_m = m + // m_E = D + // m_z = z + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); + + m_B = rateMap.hasKey("K") ? rateMap["K"].asDouble() : 0.0; + m_C = rateMap.hasKey("B") ? -rateMap["B"].asDouble() : 0.0; + m_D = rateMap.hasKey("C") ? rateMap["C"].asDouble() : 0.0; + m_m = rateMap.hasKey("m") ? rateMap["m"].asDouble() : 1.0; + m_E = rateMap.hasKey("D") ? rateMap["D"].asDouble() : 0.0; + m_z = rateMap.hasKey("z") ? rateMap["z"].asDouble() : 1.0; + m_scaling = 1.0; + + if (m_m <= 0.0 || m_z <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The Starikovskiy exponents 'm' and 'z' must be positive."); + } + } + else if (m_vibration_model == "castela") { + // Castela model: + // + // Original relaxation time: + // + // tau_k = p0 / p_k + // * exp[a_k * (T^(-1/3) - b_k) - 18.42] + // + // Equivalent bimolecular rate coefficient: + // + // k_k(T) = R T / p0 + // * exp[18.42 + a_k b_k - a_k T^(-1/3)] + // + // Internal mapping: + // + // A = R / p0 + // b = 1 + // B = 18.42 + a_k b_k + // C = -a_k + // D = 0 + // E = 0 + // scaling = 1 + + requireKey(rateMap, "a", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireKey(rateMap, "b", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "K", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + + if (rateMap.hasKey(m_reference_pressure_str)) { + m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); + } else { + m_referencePressure = OneAtm; + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Castela reference-pressure must be positive."); + } + + configureBaseFromInternalA( + node, rate_units, GasConstant / m_referencePressure, 1.0); + + m_B = 18.42 + m_castela_a * m_castela_b; + m_C = -m_castela_a; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } +} + + +void VibrationalRelaxationRate::getParameters(AnyMap& node) const +{ + if (!valid()) { + return; + } + + if (allowNegativePreExponentialFactor()) { + node["negative-A"] = true; + } + + node["vibration_model"] = m_vibration_model; + + AnyMap rateNode; + + auto storePreExponentialFactor = [&](AnyMap& target, double A) { + if (conversionUnits().factor() != 0.0) { + target[m_A_str].setQuantity(A, conversionUnits()); + } else { + target[m_A_str] = A; + target["__unconvertible__"] = true; + } + }; + + if (m_vibration_model == "constant") { + const double tol = 1e-12; + + if (std::abs(m_b) > tol + || std::abs(m_B) > tol + || std::abs(m_C) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'constant': the internal " + "parameters contain temperature-dependent terms."); + } + + storePreExponentialFactor(rateNode, m_scaling * m_A); + } + else if (m_vibration_model == "multi-state-resolved") { + storePreExponentialFactor(rateNode, m_A); + + rateNode[m_b_str] = m_b; + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; + } + else if (m_vibration_model == "starikovskiy") { + storePreExponentialFactor(rateNode, m_A); + + rateNode["n"] = m_b; + rateNode["K"] = m_B; + rateNode["B"] = -m_C; + rateNode["C"] = m_D; + rateNode["m"] = m_m; + rateNode["D"] = m_E; + rateNode["z"] = m_z; + } + else if (m_vibration_model == "castela") { + const double tol = 1e-12; + + if (std::abs(m_b - 1.0) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol + || std::abs(m_scaling - 1.0) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': the internal " + "parameters are not consistent with the Castela form. " + "Expected b = 1, D = 0, E = 0, and scaling = 1."); + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': " + "reference-pressure must be positive."); + } + + rateNode["a"] = m_castela_a; + rateNode["b"] = m_castela_b; + rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); + } + else { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } + + rateNode.setFlowStyle(); + node["rate-constant"] = std::move(rateNode); +} + + +double VibrationalRelaxationRate::ddTScaledFromStruct( + const DetailedVibData& shared_data) const +{ + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_b * invT + - (m_C / 3.0) * invT13 * invT + - m_m * m_D * std::pow(invT, m_m) * invT + - m_z * m_E * std::pow(invT, m_z) * invT; +} + + +void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + if (rxn.reversible) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Vibrational relaxation rates do not support reversible " + "reactions."); + } + + const string family = inferRelaxingFamily(rxn); + + if (m_vibration_model == "constant") { + validateSimpleRelaxationToGroundState(rxn, "constant"); + } else if (m_vibration_model == "castela") { + validateCastelaReaction(rxn); + } else if (m_vibration_model == "starikovskiy") { + validateSimpleRelaxationToGroundState(rxn, "starikovskiy"); + } else if (m_vibration_model == "multi-state-resolved") { + validateDetailedRelaxationReaction(rxn); + } else { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Unrecognized vibration_model '{}'.", m_vibration_model); + } + + registerVibrationalModelConsistency( + kin, family, m_vibration_model, rxn.input); +} + +} // namespace Cantera \ No newline at end of file diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 7fa38c916e5..1961fa2a964 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -21,7 +21,7 @@ typedef Eigen::SparseMatrix SparseMat; EEDFTwoTermApproximation::EEDFTwoTermApproximation(PlasmaPhase* s) { - // store a pointer to s. + // Store the PlasmaPhase context used by the solver (pointer to s). m_phase = s; m_first_call = true; m_has_EEDF = false; @@ -43,6 +43,96 @@ void EEDFTwoTermApproximation::setLinearGrid(double& kTe_max, size_t& ncell) setGridCache(); } +void EEDFTwoTermApproximation::setQuadraticGrid(double& kTe_max, size_t& ncell) +{ + m_points = ncell; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + double n = static_cast(m_points); + + for (size_t j = 0; j <= m_points; j++) { + double x = static_cast(j); + m_gridEdge[j] = kTe_max * x * (x + 1.0) / (n * (n + 1.0)); + } + + for (size_t j = 0; j < m_points; j++) { + m_gridCenter[j] = 0.5 * (m_gridEdge[j] + m_gridEdge[j + 1]); + } + + setGridCache(); +} + +void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell, double ratio) +{ + m_points = ncell; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + double denominator = std::pow(ratio, m_points) - 1.0; + + if (std::abs(denominator) < 1e-14) { + throw CanteraError("EEDFTwoTermApproximation::setGeometricGrid", + "Invalid geometric-grid ratio."); + } + + for (size_t j = 0; j < m_points; j++) { + m_gridEdge[j] = kTe_max * (std::pow(ratio, j) - 1.0) / denominator; + + m_gridCenter[j] = kTe_max + * (std::pow(ratio, j + 0.5) - 1.0) + / denominator; + } + + m_gridEdge[m_points] = kTe_max; + + setGridCache(); +} + +void EEDFTwoTermApproximation::setCustomGrid(span levels) +{ + if (levels.size() < 2) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid must contain at least two edge points."); + } + + m_points = levels.size() - 1; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + for (size_t j = 0; j < m_points + 1; j++) { + if (!std::isfinite(levels[j])) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid contains a non-finite value."); + } + if (levels[j] < 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid values must be non-negative."); + } + if (j > 0 && levels[j] <= levels[j - 1]) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid values must be strictly increasing."); + } + + m_gridEdge[j] = levels[j]; + } + + for (size_t j = 0; j < m_points; j++) { + m_gridCenter[j] = 0.5 * (m_gridEdge[j] + m_gridEdge[j + 1]); + } + + setGridCache(); +} + int EEDFTwoTermApproximation::calculateDistributionFunction() { if (m_first_call) { @@ -51,33 +141,133 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() } updateMoleFractions(); + checkSpeciesNoCrossSection(); updateCrossSections(); - if (!m_has_EEDF) { - if (m_firstguess == "maxwell") { - for (size_t j = 0; j < m_points; j++) { - m_f0(j) = 2.0 * pow(1.0 / Pi, 0.5) * pow(m_init_kTe, -3. / 2.) * - exp(-m_gridCenter[j] / m_init_kTe); + const double EN = m_phase->reducedElectricField(); + + // Helper used to impose a Maxwellian EEDF. + auto setMaxwellian = [&](double kTe_eV) { + if (!std::isfinite(kTe_eV) || kTe_eV <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Invalid kTe value for Maxwellian first guess."); + } + + const double fFloor = 1e-300; + + for (size_t j = 0; j < m_points; j++) { + double arg = -m_gridCenter[j] / kTe_eV; + + if (arg < -700.0) { + m_f0(j) = fFloor; + } else { + m_f0(j) = std::max(fFloor, + 2.0 * std::sqrt(1.0 / Pi) * std::pow(kTe_eV, -1.5) + * std::exp(arg)); } - } else { + } + + double fnorm = norm(m_f0, m_gridCenter); + + if (!std::isfinite(fnorm) || fnorm <= 0.0) { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - " unknown EEDF first guess"); + "Invalid norm for Maxwellian initialization."); } - } - converge(m_f0); + m_f0 /= fnorm; + }; + + // At very low reduced electric field, force a Maxwellian at the gas + // temperature and skip the Boltzmann convergence to avoid wrong EEDF convergence and numerical instabilities + // when integrating over time. + if (EN <= EN_min) { + setMaxwellian(Boltzmann * m_phase->temperature() / ElectronCharge); + } else { + // At non-zero reduced electric field, use a hot Maxwellian first guess + // only when no previously converged EEDF is available. + if (!m_has_EEDF) { + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + } + + converge(m_f0); + + // Grid adaptation based on EEDF tail decay. If enabled, this will iteratively adjust the grid + // until the EEDF tail decays within the specified bounds given in the YAML file. + // @todo Implement a more robust version which also varies the number of grid points if necessary. + // The current version only adjusts the maximum energy of the grid. + + if (m_adaptGrid) { + const double fFloor = 1e-300; + + for (size_t n = 0; n < m_maxGridAdaptIterations; n++) { + double fLeft = std::max(std::abs(m_f0(0)), fFloor); + double fRight = std::max(std::abs(m_f0(m_points - 1)), fFloor); + + double decades = std::log10(fLeft) - std::log10(fRight); + + if (!std::isfinite(decades)) { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Non-finite EEDF decay detected during grid adaptation."); + } + + if (decades < m_minEedfDecay) { + // The right boundary is too low: the tail has not decayed enough. + double newMaxEnergy = m_kTeMax * (1.0 + m_gridUpdateFactor); + + updateGrid(newMaxEnergy); + + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + + updateCrossSections(); + converge(m_f0); - // write the EEDF at grid edges + } else if (decades > m_maxEedfDecay) { + // The right boundary is unnecessarily high. + double newMaxEnergy = m_kTeMax / (1.0 + m_gridUpdateFactor); + + updateGrid(newMaxEnergy); + + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + + updateCrossSections(); + converge(m_f0); + + } else { + break; + } + } + } + } + + // Write the EEDF at grid edges. vector f(m_f0.data(), m_f0.data() + m_f0.rows() * m_f0.cols()); - vector x(m_gridCenter.data(), m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); + vector x(m_gridCenter.data(), + m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); + for (size_t i = 0; i < m_points + 1; i++) { m_f0_edge[i] = linearInterp(m_gridEdge[i], x, f); } m_has_EEDF = true; - // update electron mobility - m_electronMobility = electronMobility(m_f0); + // Update electron mobility. + m_electronMobility = computeElectronMobility(m_f0); + return 0; } @@ -115,7 +305,7 @@ void EEDFTwoTermApproximation::converge(Eigen::VectorXd& f0) if (err1 < m_rtol) { break; } else if (n == m_maxn - 1) { - throw CanteraError("WeaklyIonizedGas::converge", "Convergence failed"); + throw CanteraError("EEDFTwoTermApproximation::converge", "Convergence failed"); } } } @@ -132,7 +322,7 @@ Eigen::VectorXd EEDFTwoTermApproximation::iterate(const Eigen::VectorXd& f0, dou for (size_t k : m_phase->kInelastic()) { SparseMat Q_k = matrix_Q(g, k); SparseMat P_k = matrix_P(g, k); - PQ += (matrix_Q(g, k) - matrix_P(g, k)) * m_X_targets[m_klocTargets[k]]; + PQ += (Q_k - P_k) * m_X_targets[m_klocTargets[k]]; } SparseMat A = matrix_A(f0); @@ -278,7 +468,7 @@ SparseMat EEDFTwoTermApproximation::matrix_A(const Eigen::VectorXd& f0) double alpha; double E = m_phase->electricField(); if (m_growth == "spatial") { - double mu = electronMobility(f0); + double mu = computeElectronMobility(f0); double D = electronDiffusivity(f0); alpha = (mu * E - sqrt(pow(mu * E, 2) - 4 * D * nu * nDensity)) / 2.0 / D / nDensity; } else { @@ -342,7 +532,7 @@ SparseMat EEDFTwoTermApproximation::matrix_A(const Eigen::VectorXd& f0) SparseMat A(m_points, m_points); A.setFromTriplets(tripletList.begin(), tripletList.end()); - //plus G + // plus G SparseMat G(m_points, m_points); if (m_growth == "temporal") { for (size_t i = 0; i < m_points; i++) { @@ -394,7 +584,7 @@ double EEDFTwoTermApproximation::electronDiffusivity(const Eigen::VectorXd& f0) return 1./3. * m_gamma * simpson(f, x) / nDensity; } -double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) +double EEDFTwoTermApproximation::computeElectronMobility(const Eigen::VectorXd& f0) { double nu = netProductionFrequency(f0); vector y(m_points + 1, 0.0); @@ -407,7 +597,9 @@ double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) } } double nDensity = m_phase->molarDensity() * Avogadro; - return -1./3. * m_gamma * simpson(asVectorXd(y), asVectorXd(m_gridEdge)) / nDensity; + auto f = ConstMappedVector(y.data(), y.size()); + auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size()); + return -1./3. * m_gamma * simpson(f, x) / nDensity; } void EEDFTwoTermApproximation::initSpeciesIndexCrossSections() @@ -464,6 +656,7 @@ void EEDFTwoTermApproximation::updateCrossSections() } // Update the species mole fractions used for EEDF computation +// Renormalize over species with electron-collision cross-section data. void EEDFTwoTermApproximation::updateMoleFractions() { double tmp_sum = 0.0; @@ -478,38 +671,144 @@ void EEDFTwoTermApproximation::updateMoleFractions() } } +// The former implementation counted the effective cros section as an elastic contribution +// thus double counting the inelastic contributions to the total cross section. +// This implemetation avoids this drawback. void EEDFTwoTermApproximation::calculateTotalCrossSection() { m_totalCrossSectionCenter.assign(m_points, 0.0); m_totalCrossSectionEdge.assign(m_points + 1, 0.0); + + std::vector has_effective_for_target(m_phase->nSpecies(), false); + + // build this list first to avoid scanning through all collisions for each target twice in the loop. + for (size_t ke : m_phase->kElastic()) { + if (m_phase->collisionRate(ke)->kind() == "effective") { + has_effective_for_target[m_phase->targetIndex(ke)] = true; + } + } + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + const std::string& kind = m_phase->collisionRate(k)->kind(); + const size_t target = m_phase->targetIndex(k); + + // If this target has an effective cross section, then that effective + // cross section already contains elastic + inelastic contributions. + // So only the effective cross section contributes to sigma_total. + if (has_effective_for_target[target] && kind != "effective") { + continue; + } + auto x = m_phase->collisionRate(k)->energyLevels(); auto y = m_phase->collisionRate(k)->crossSections(); + const double X_target = m_X_targets[m_klocTargets[k]]; + for (size_t i = 0; i < m_points; i++) { - m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * + m_totalCrossSectionCenter[i] += X_target * linearInterp(m_gridCenter[i], x, y); } + for (size_t i = 0; i < m_points + 1; i++) { - m_totalCrossSectionEdge[i] += m_X_targets[m_klocTargets[k]] * + m_totalCrossSectionEdge[i] += X_target * linearInterp(m_gridEdge[i], x, y); } } } +// new implementation of the function proposed by nicolas void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() { m_sigmaElastic.clear(); m_sigmaElastic.resize(m_points, 0.0); for (size_t k : m_phase->kElastic()) { + + writelog("Enter with an elastic cs\n"); auto x = m_phase->collisionRate(k)->energyLevels(); auto y = m_phase->collisionRate(k)->crossSections(); + vector y_elastic(y.data(), y.data() + y.size()); + // Note: // moleFraction(m_kTargets[k]) <=> m_X_targets[m_klocTargets[k]] double mass_ratio = ElectronMass / (m_phase->molecularWeight(m_kTargets[k]) / Avogadro); + + if (m_phase->collisionRate(k)->kind()=="effective") { + + writelog("Enter with an elastic cs that is actually an effective\n"); + + for (size_t ki : m_phase->kInelastic()) + { + if(m_phase->targetIndex(ki) == m_phase->targetIndex(k)){ + writelog("loop over inelastic processes: process {}\n", ki); + auto xi = m_phase->collisionRate(ki)->energyLevels(); + auto yi = m_phase->collisionRate(ki)->crossSections(); + for (size_t i = 0; i < x.size(); i++) + { + y_elastic[i] -= linearInterpCrossSectionZeroOutside(x[i], xi, yi); + } + } + } + // check that the reconstructed elastic cross section is non-negative. + for (size_t i = 0; i < y_elastic.size(); i++) { + if (y_elastic[i] < 0.0) { + if (y_elastic[i] > -1e-30) { + y_elastic[i] = 0.0; + } else { + const std::string& effectiveKind = m_phase->collisionRate(k)->kind(); + std::string effectiveTarget = m_phase->speciesName(m_phase->targetIndex(k)); + + writelog("Warning: reconstructed elastic cross section is negative " + "for collision {} kind {} target {} at energy {}: {}\n", + k, effectiveKind, effectiveTarget, x[i], y_elastic[i]); + y_elastic[i] = 0.0; + } + } + } + } + for (size_t i = 0; i < m_points; i++) { m_sigmaElastic[i] += 2.0 * mass_ratio * m_X_targets[m_klocTargets[k]] * - linearInterp(m_gridEdge[i], x, y); + linearInterp(m_gridEdge[i], x, y_elastic); + } + } +} + +double EEDFTwoTermApproximation::linearInterpBounded(double x, + span xpts, + span fpts, + double below_value, + double above_value) +{ + AssertThrowMsg(!xpts.empty(), "linearInterpBounded", "x data empty"); + AssertThrowMsg(!fpts.empty(), "linearInterpBounded", "f(x) data empty"); + AssertThrowMsg(xpts.size() == fpts.size(), "linearInterpBounded", + "len(xpts) = {}, len(fpts) = {}", xpts.size(), fpts.size()); + + if (x < xpts.front()) { + return below_value; + } + + if (x > xpts.back()) { + return above_value; + } + + return linearInterp(x, xpts, fpts); +} + +double EEDFTwoTermApproximation::linearInterpCrossSectionZeroOutside(double x, + span xpts, + span fpts) +{ + return linearInterpBounded(x, xpts, fpts, 0.0, 0.0); +} + +void EEDFTwoTermApproximation::checkSpeciesNoCrossSection() +{ + // warn that a specific species needs cross-section data. + for (size_t k : m_kOthers) { + if (m_phase->moleFraction(k) > m_moleFractionThreshold) { + writelog("EEDFTwoTermApproximation:checkSpeciesNoCrossSection\n"); + writelog("Warning:The mole fraction of species {} is more than 0.01 (X = {:.3g}) but it has no cross-section data\n", m_phase->speciesName(k), m_phase->moleFraction(k)); } } } @@ -596,4 +895,111 @@ double EEDFTwoTermApproximation::norm(const Eigen::VectorXd& f, const Eigen::Vec return numericalQuadrature(m_quadratureMethod, p, grid); } +void EEDFTwoTermApproximation::setGridType(const string& gridType) +{ + if (gridType != "Linear" && + gridType != "Quadratic" && + gridType != "Geometric") { + throw CanteraError("EEDFTwoTermApproximation::setGridType", + "Unknown energy grid type '{}'. Expected Linear, Quadratic or Geometric.", + gridType); + } + + m_gridType = gridType; +} + +void EEDFTwoTermApproximation::setInitialGridParameters(double initialMaxEnergy, + size_t nGridCells) +{ + if (!std::isfinite(initialMaxEnergy) || initialMaxEnergy <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setInitialGridParameters", + "initialMaxEnergy must be finite and greater than zero."); + } + + if (nGridCells == 0) { + throw CanteraError("EEDFTwoTermApproximation::setInitialGridParameters", + "nGridCells must be greater than zero."); + } + + m_kTeMax = initialMaxEnergy; + m_initialGridCells = nGridCells; +} + +void EEDFTwoTermApproximation::enableGridAdaptation(bool enabled) +{ + m_adaptGrid = enabled; +} + +void EEDFTwoTermApproximation::setGridAdaptationParameters(bool enabled, + double minDecayDecades, + double maxDecayDecades, + double updateFactor, + size_t maxIterations) +{ + if (!std::isfinite(minDecayDecades) || !std::isfinite(maxDecayDecades) || + minDecayDecades <= 0.0 || maxDecayDecades <= minDecayDecades) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "Require 0 < min_decay_decades < max_decay_decades."); + } + + if (!std::isfinite(updateFactor) || updateFactor <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "update_factor must be finite and greater than zero."); + } + + if (maxIterations == 0) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "max_iterations must be greater than zero."); + } + + m_adaptGrid = enabled; + m_minEedfDecay = minDecayDecades; + m_maxEedfDecay = maxDecayDecades; + m_gridUpdateFactor = updateFactor; + m_maxGridAdaptIterations = maxIterations; +} + +void EEDFTwoTermApproximation::updateGrid(double maxEnergy) +{ + if (!std::isfinite(maxEnergy) || maxEnergy <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::updateGrid", + "Maximum grid energy must be finite and greater than zero."); + } + + m_kTeMax = maxEnergy; + + if (m_gridType == "Linear") { + setLinearGrid(m_kTeMax, m_initialGridCells); + } else if (m_gridType == "Quadratic") { + setQuadraticGrid(m_kTeMax, m_initialGridCells); + } else if (m_gridType == "Geometric") { + setGeometricGrid(m_kTeMax, m_initialGridCells); + } else { + throw CanteraError("EEDFTwoTermApproximation::updateGrid", + "Unknown energy grid type '{}'.", m_gridType); + } + + // put back the flag at false because since the grid has changed the EEDF must be computed again. + m_has_EEDF = false; +} + +double EEDFTwoTermApproximation::getElectronMobility() const +{ + if (!m_has_EEDF || !std::isfinite(m_electronMobility)) { + throw CanteraError("EEDFTwoTermApproximation::getElectronMobility", + "Electron mobility is not available before a valid EEDF has been computed."); + } + return m_electronMobility; +} + +// Set the reduced electric field threshold below which the EEDF is forced to be Maxwellian at the gas temperature. +// The input to this function is expected to be in Townsend. +void EEDFTwoTermApproximation::setReducedFieldThresholdBeforeMaxwellianTd(double threshold){ + if (!std::isfinite(threshold) || threshold < 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setReducedFieldThresholdBeforeMaxwellianTd", + "Reduced field threshold must be finite and non-negative."); + } + EN_min = threshold*1e-21; +} + } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 26f83aeb469..28b9ee95d84 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -19,25 +19,56 @@ namespace Cantera { namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); + + bool isVibrationalReservoirName(const string& name, string& baseName) + { + const string suffix = "(v)"; + + if (name.size() <= suffix.size()) { + return false; + } + + if (name.compare(name.size() - suffix.size(), suffix.size(), suffix) != 0) { + return false; + } + + baseName = name.substr(0, name.size() - suffix.size()); + return !baseName.empty(); + } } PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { - initThermoFile(inputFile, id_); - - // initial electron temperature + // Initialize electron temperature before setParameters() can trigger EEDF updates. m_electronTemp = temperature(); - // Initialize the Boltzmann Solver + // The EEDF solver must exist before initThermoFile(), because setParameters() + // may add electron collisions and addCollision() updates the solver cache. m_eedfSolver = make_unique(this); - // Set Energy Grid (Hardcoded Defaults for Now) - double kTe_max = 60; - size_t nGridCells = 301; - m_nPoints = nGridCells + 1; + // Safe default grid used before / unless the input file overrides it. + double kTe_max = 100.0; + size_t nGridCells = 1001; + m_eedfSolver->setInitialGridParameters(kTe_max, nGridCells); m_eedfSolver->setLinearGrid(kTe_max, nGridCells); - m_electronEnergyLevels = asVectorXd(m_eedfSolver->getGridEdge()); - m_electronEnergyDist.setZero(m_nPoints); + + auto levels = m_eedfSolver->getGridEdge(); + m_nPoints = levels.size(); + m_electronEnergyLevels = Eigen::Map( + levels.data(), m_nPoints); + m_electronEnergyDist.resize(m_nPoints); + m_electronEnergyDist.setZero(); + + if (!inputFile.empty()) { + initThermoFile(inputFile, id_); + } + + // If no EEDF was supplied by input, initialize a valid default isotropic EEDF. + if (m_distributionType == "isotropic" && + static_cast(m_electronEnergyDist.size()) == m_nPoints && + m_electronEnergyDist.sum() == 0.0) { + updateElectronEnergyDistribution(); + } } PlasmaPhase::~PlasmaPhase() @@ -60,10 +91,28 @@ void PlasmaPhase::updateElectronEnergyDistribution() } else if (m_distributionType == "isotropic") { setIsotropicElectronEnergyDistribution(); } else if (m_distributionType == "Boltzmann-two-term") { + // The two-term Boltzmann solver may adapt the electron energy grid, so both + // the grid and the distribution are copied back before invalidating dependent + // cross-section data. auto ierr = m_eedfSolver->calculateDistributionFunction(); if (ierr == 0) { + auto levels = m_eedfSolver->getGridEdge(); auto y = m_eedfSolver->getEEDFEdge(); - m_electronEnergyDist = Eigen::Map(y.data(), m_nPoints); + + if (levels.size() != y.size()) { + throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", + "Inconsistent EEDF solver output: grid edge size and EEDF edge size differ."); + } + + m_nPoints = levels.size(); + + m_electronEnergyLevels = Eigen::Map( + levels.data(), m_nPoints); + + m_electronEnergyDist = Eigen::Map( + y.data(), m_nPoints); + + electronEnergyLevelChanged(); } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Call to calculateDistributionFunction failed."); @@ -78,6 +127,9 @@ void PlasmaPhase::updateElectronEnergyDistribution() if (validEEDF) { updateElectronTemperatureFromEnergyDist(); } else { + // Keep the previous electron temperature if the solver did not return a usable + // distribution. This avoids replacing Te with a value derived from invalid data. + // The user will be warned of this behaviour. writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); } } else { @@ -88,16 +140,22 @@ void PlasmaPhase::updateElectronEnergyDistribution() electronEnergyDistributionChanged(); } -void PlasmaPhase::normalizeElectronEnergyDistribution() { - Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3./2.); - double norm = 2./3. * numericalQuadrature(m_quadratureMethod, - m_electronEnergyDist, eps32); - if (norm < 0.0) { +void PlasmaPhase::normalizeElectronEnergyDistribution() +{ + checkElectronEnergyLevels(); + checkElectronEnergyDistribution(); + + Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3.0 / 2.0); + double norm = 2.0 / 3.0 * numericalQuadrature( + m_quadratureMethod, m_electronEnergyDist, eps32); + + if (!std::isfinite(norm) || norm <= 0.0) { throw CanteraError("PlasmaPhase::normalizeElectronEnergyDistribution", - "The norm is negative. This might be caused by bad " - "electron energy distribution"); + "Electron energy distribution has invalid norm: {}.", norm); } + m_electronEnergyDist /= norm; + checkElectronEnergyDistribution(); } void PlasmaPhase::setElectronEnergyDistributionType(const string& type) @@ -127,7 +185,13 @@ void PlasmaPhase::setIsotropicElectronEnergyDistribution() checkElectronEnergyDistribution(); } -void PlasmaPhase::setElectronTemperature(const double Te) { +void PlasmaPhase::setElectronTemperature(const double Te) +{ + if (!std::isfinite(Te) || Te <= 0.0) { + throw CanteraError("PlasmaPhase::setElectronTemperature", + "Electron temperature must be finite and positive."); + } + m_electronTemp = Te; updateElectronEnergyDistribution(); } @@ -155,7 +219,13 @@ void PlasmaPhase::endEquilibrate() ThermoPhase::endEquilibrate(); } -void PlasmaPhase::setMeanElectronEnergy(double energy) { +void PlasmaPhase::setMeanElectronEnergy(double energy) +{ + if (!std::isfinite(energy) || energy <= 0.0) { + throw CanteraError("PlasmaPhase::setMeanElectronEnergy", + "Mean electron energy must be finite and positive."); + } + m_electronTemp = 2.0 / 3.0 * energy * ElectronCharge / Boltzmann; updateElectronEnergyDistribution(); } @@ -179,49 +249,102 @@ void PlasmaPhase::electronEnergyLevelChanged() m_levelNum++; // Cross sections are interpolated on the energy levels if (m_collisions.size() > 0) { + vector energyLevels(m_nPoints); + MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; for (shared_ptr collision : m_collisions) { const auto& rate = boost::polymorphic_pointer_downcast (collision->rate()); - rate->updateInterpolatedCrossSection(asSpan(m_electronEnergyLevels)); + rate->updateInterpolatedCrossSection(energyLevels); } } } void PlasmaPhase::checkElectronEnergyLevels() const { - Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) - - m_electronEnergyLevels.head(m_nPoints - 1); - if (m_electronEnergyLevels[0] < 0.0 || (h <= 0.0).any()) { + if (m_nPoints < 2 || + static_cast(m_electronEnergyLevels.size()) != m_nPoints) { throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", - "Values of electron energy levels need to be positive and " - "monotonically increasing."); + "Electron energy levels must contain at least two points and match m_nPoints."); + } + + for (size_t i = 0; i < m_nPoints; i++) { + const double eps = m_electronEnergyLevels[i]; + if (!std::isfinite(eps)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy level {} is non-finite.", i); + } + if (eps < 0.0) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy levels must be non-negative."); + } + if (i > 0 && eps <= m_electronEnergyLevels[i - 1]) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy levels must be strictly increasing."); + } } } + void PlasmaPhase::checkElectronEnergyDistribution() const { - Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) - - m_electronEnergyLevels.head(m_nPoints - 1); - if ((m_electronEnergyDist < 0.0).any()) { + if (m_nPoints < 2 || + static_cast(m_electronEnergyDist.size()) != m_nPoints || + static_cast(m_electronEnergyLevels.size()) != m_nPoints) { throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", - "Values of electron energy distribution cannot be negative."); + "Electron energy distribution and energy levels must have matching " + "sizes of at least two points."); } + + double sum = 0.0; + double maxVal = -BigNumber; + + for (size_t i = 0; i < m_nPoints; i++) { + const double f = m_electronEnergyDist[i]; + if (!std::isfinite(f)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution contains a non-finite value at index {}.", + i); + } + if (f < 0.0) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution cannot contain negative values."); + } + sum += f; + maxVal = std::max(maxVal, f); + } + + if (!(sum > 0.0) || !(maxVal > 0.0)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution must contain at least one positive value."); + } + if (m_electronEnergyDist[m_nPoints - 1] > 0.01) { warn_user("PlasmaPhase::checkElectronEnergyDistribution", - "The value of the last element of electron energy distribution exceed 0.01. " - "This indicates that the value of electron energy level is not high enough " - "to contain the isotropic distribution at mean electron energy of " - "{} eV", meanElectronEnergy()); + "The last value of the electron energy distribution exceeds 0.01. " + "The electron energy grid may be too short for a mean electron energy " + "of {} eV.", meanElectronEnergy()); } } void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, span dist) { + if (levels.size() != dist.size()) { + throw CanteraError("PlasmaPhase::setDiscretizedElectronEnergyDist", + "Energy levels and electron energy distribution must have the same size. " + "Got {} levels and {} distribution values.", levels.size(), dist.size()); + } + + if (levels.size() < 2) { + throw CanteraError("PlasmaPhase::setDiscretizedElectronEnergyDist", + "A discretized electron energy distribution requires at least two points."); + } m_distributionType = "discretized"; m_nPoints = levels.size(); - m_electronEnergyLevels = asVectorXd(levels); - m_electronEnergyDist = asVectorXd(dist); + m_electronEnergyLevels = + Eigen::Map(levels.data(), m_nPoints); + m_electronEnergyDist = + Eigen::Map(dist.data(), m_nPoints); checkElectronEnergyLevels(); if (m_do_normalizeElectronEnergyDist) { normalizeElectronEnergyDistribution(); @@ -245,15 +368,22 @@ void PlasmaPhase::updateElectronTemperatureFromEnergyDist() "trapezoidal", m_electronEnergyDist, eps52); } - if (epsilon_m < 0.0) { + if (!std::isfinite(epsilon_m) || epsilon_m <= 0.0) { throw CanteraError("PlasmaPhase::updateElectronTemperatureFromEnergyDist", - "The electron energy distribution produces negative electron temperature."); + "The electron energy distribution produces an invalid mean electron energy: {}.", + epsilon_m); } m_electronTemp = 2.0 / 3.0 * epsilon_m * ElectronCharge / Boltzmann; } -void PlasmaPhase::setIsotropicShapeFactor(double x) { +void PlasmaPhase::setIsotropicShapeFactor(double x) +{ + if (!std::isfinite(x) || x <= 0.0) { + throw CanteraError("PlasmaPhase::setIsotropicShapeFactor", + "The isotropic shape factor must be finite and positive."); + } + m_isotropicShapeFactor = x; updateElectronEnergyDistribution(); } @@ -318,32 +448,242 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) auto levels = eedf["energy-levels"].asVector(); auto distribution = eedf["distribution"].asVector(levels.size()); setDiscretizedElectronEnergyDist(levels, distribution); + } else if (m_distributionType == "Boltzmann-two-term") { + m_eedfSolver = make_unique(this); + + if (eedf.hasKey("energy-levels")) { + // Mode A: user-provided grid edges. + // In this mode, the grid is considered explicit and fixed by default. + auto levels = eedf["energy-levels"].asVector(); + + m_eedfSolver->setCustomGrid(levels); + m_eedfSolver->enableGridAdaptation(false); + + m_nPoints = levels.size(); + + } else { + // Mode B: generated initial grid. + if (!eedf.hasKey("initial_max_energy_level")) { + throw CanteraError("PlasmaPhase::setParameters", + "Boltzmann-two-term requires either 'energy-levels' or " + "'initial_max_energy_level'."); + } + + if (!eedf.hasKey("initial_number_of_energy_grid_cells")) { + throw CanteraError("PlasmaPhase::setParameters", + "Boltzmann-two-term requires either 'energy-levels' or " + "'initial_number_of_energy_grid_cells'."); + } + + double initialMaxEnergy = eedf["initial_max_energy_level"].asDouble(); + size_t nGridCells = static_cast( + eedf["initial_number_of_energy_grid_cells"].asInt()); + + if (!std::isfinite(initialMaxEnergy) || initialMaxEnergy <= 0.0) { + throw CanteraError("PlasmaPhase::setParameters", + "initial_max_energy_level must be finite and greater than zero."); + } + + if (nGridCells == 0) { + throw CanteraError("PlasmaPhase::setParameters", + "initial_number_of_energy_grid_cells must be greater than zero."); + } + + string energyLevelsDistribution = "Linear"; + if (eedf.hasKey("energy_levels_distribution")) { + energyLevelsDistribution = + eedf["energy_levels_distribution"].asString(); + } else { + writelog("No energy_levels_distribution key found in the input file. " + "Defaulting to linear grid.\n"); + } + + m_eedfSolver->setGridType(energyLevelsDistribution); + m_eedfSolver->setInitialGridParameters(initialMaxEnergy, nGridCells); + + if (energyLevelsDistribution == "Linear") { + m_eedfSolver->setLinearGrid(initialMaxEnergy, nGridCells); + } else if (energyLevelsDistribution == "Quadratic") { + m_eedfSolver->setQuadraticGrid(initialMaxEnergy, nGridCells); + } else if (energyLevelsDistribution == "Geometric") { + if (eedf.hasKey("geometric_grid_ratio")) { + double ratio = eedf["geometric_grid_ratio"].asDouble(); + if (!std::isfinite(ratio) || ratio <= 1.0) { + throw CanteraError("PlasmaPhase::setParameters", + "geometric_grid_ratio must be finite and greater than 1.0."); + } + m_eedfSolver->setGeometricGrid(initialMaxEnergy, nGridCells, ratio); + } else { + writelog("No geometric_grid_ratio key found for geometric grid in the input file. " + "Defaulting to a geometric ratio of 1.01.\n"); + m_eedfSolver->setGeometricGrid(initialMaxEnergy, nGridCells); + } + } else { + throw CanteraError("PlasmaPhase::setParameters", + "energy_levels_distribution should be Linear, Quadratic or Geometric."); + } + + if (eedf.hasKey("reduced_field_threshold_before_maxwellian_Td")){ + double maxwellian_threshold = eedf["reduced_field_threshold_before_maxwellian_Td"].asDouble(); + if (!std::isfinite(maxwellian_threshold) || maxwellian_threshold < 0.0) { + throw CanteraError("PlasmaPhase::setParameters", + "reduced_field_threshold_before_maxwellian_Td must be finite and non-negative."); + } + m_eedfSolver->setReducedFieldThresholdBeforeMaxwellianTd(maxwellian_threshold); // The input to this function is expected to be in Townsend. + } + + if (eedf.hasKey("energy_grid_adaptation")) { + const AnyMap adapt = eedf["energy_grid_adaptation"].as(); + + bool enabled = true; + if (adapt.hasKey("enabled")) { + enabled = adapt["enabled"].asBool(); + } + + double minDecayDecades = 8.0; + if (adapt.hasKey("min_decay_decades")) { + minDecayDecades = adapt["min_decay_decades"].asDouble(); + } + + double maxDecayDecades = 14.0; + if (adapt.hasKey("max_decay_decades")) { + maxDecayDecades = adapt["max_decay_decades"].asDouble(); + } + + double updateFactor = 0.25; + if (adapt.hasKey("update_factor")) { + updateFactor = adapt["update_factor"].asDouble(); + } + + size_t maxIterations = 5; + if (adapt.hasKey("max_iterations")) { + maxIterations = static_cast( + adapt["max_iterations"].asInt()); + } + + m_eedfSolver->setGridAdaptationParameters( + enabled, minDecayDecades, maxDecayDecades, + updateFactor, maxIterations); + } else { + m_eedfSolver->enableGridAdaptation(false); + } + + // In PlasmaPhase, m_nPoints is the number of grid edges. + m_nPoints = nGridCells + 1; + } + + m_electronEnergyLevels = Eigen::Map( + m_eedfSolver->getGridEdge().data(), m_eedfSolver->getGridEdge().size()); + + m_nPoints = m_eedfSolver->getGridEdge().size(); + + m_electronEnergyDist.resize(m_nPoints); + m_electronEnergyDist.setZero(); + + checkElectronEnergyLevels(); + electronEnergyLevelChanged(); } + + } if (rootNode.hasKey("electron-collisions")) { for (const auto& item : rootNode["electron-collisions"].asVector()) { - auto rate = make_shared(item); - Composition reactants, products; - reactants[item["target"].asString()] = 1; - reactants[electronSpeciesName()] = 1; - if (item.hasKey("product")) { - products[item["product"].asString()] = 1; - } else { - products[item["target"].asString()] = 1; + if (item.hasKey("name")) { + string name = item["name"].asString(); + if (name.empty()) { + throw InputFileError("setParameters", rootNode, "Empty electron-collision name."); + } + if (m_electronCollisionDefinitions.count(name)) { + throw InputFileError("setParameters", rootNode, "Duplicate electron-collision name '{}'.", name); + } + m_electronCollisionDefinitions[name] = item; } - products[electronSpeciesName()] = 1; - if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; - } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; + } + } + + // in the case of a wrong combination of the two entry formats (should the user mix new and old formats somehow), ensure that each collision is only loaded once + if (rootNode.hasKey("reactions")) { + for (const auto& rxnNode : rootNode["reactions"].asVector()) { + if (rxnNode.hasKey("type") + && rxnNode["type"].asString() == "electron-collision-plasma" + && rxnNode.hasKey("collision")) { + string name = rxnNode["collision"].asString(); + + if (!m_electronCollisionDefinitions.count(name)) { + throw InputFileError("setParameters", rootNode, + "Reaction references unknown electron collision '{}'.", name); + } + m_referencedElectronCollisions.insert(name); } - auto R = make_shared(reactants, products, rate); - addCollision(R); } } + + if (rootNode.hasKey("electron-collisions")) { + for (const auto& item : rootNode["electron-collisions"].asVector()) { + if (item.hasKey("name")) { + string name = item["name"].asString(); + + if (m_referencedElectronCollisions.count(name)) { + continue; + } + + } + bool hasName = item.hasKey("name"); + + if (hasName) { + string name = item["name"].asString(); + + // New YAML format: if the collision is referenced by a reaction, do not create a synthetic reaction here. + if (m_referencedElectronCollisions.count(name)) { + continue; + } + } + + // Old YAML format for anonymous collisions or unreferenced named collisions: they remain as standalone EEDF collisions. + addStandaloneElectronCollision(item); + + } + } +} + + +void PlasmaPhase::addStandaloneElectronCollision(const AnyMap& item) +{ + auto rate = make_shared(item); + + if (!rate->hasCrossSectionData()) { + throw InputFileError("addStandaloneElectronCollision", item, + "Cross-sections are required in the reaction if you are using the deprecated input format!"); + } + + Composition reactants; + Composition products; + + string target = item["target"].asString(); + + reactants[target] = 1; + reactants[electronSpeciesName()] = 1; + + if (item.hasKey("product")) { + products[item["product"].asString()] = 1; + } else { + products[target] = 1; + } + + products[electronSpeciesName()] = 1; + + if (rate->kind() == "ionization") { + products[electronSpeciesName()] += 1; + } else if (rate->kind() == "attachment") { + products[electronSpeciesName()] -= 1; + } + + auto R = make_shared(reactants, products, rate); + addCollision(R); } + bool PlasmaPhase::addSpecies(shared_ptr spec) { bool added = IdealGasPhase::addSpecies(spec); @@ -361,6 +701,11 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) "Only one electron species is allowed.", spec->name); } } + + if (added) { + m_vibrationalReservoirSpeciesNeedUpdate = true; + } + return added; } @@ -422,6 +767,7 @@ void PlasmaPhase::setCollisions() void PlasmaPhase::addCollision(shared_ptr collision) { + size_t i = nCollisions(); // setup callback to signal updating the cross-section-related @@ -447,6 +793,63 @@ void PlasmaPhase::addCollision(shared_ptr collision) " collision with equation '{}'", collision->equation()); } + // management of the new data file format + + auto ratePtr = std::dynamic_pointer_cast(collision->rate()); + + if (!ratePtr) { + throw CanteraError("addCollision", "The rate is not initialised"); + } + + if (!ratePtr->collisionName().empty() && !ratePtr->hasCrossSectionData()) { + auto it = m_electronCollisionDefinitions.find(ratePtr->collisionName()); + if (it == m_electronCollisionDefinitions.end()) { + throw CanteraError("addCollision", "Unknown electron collision '{}'.", ratePtr->collisionName()); + } + + ratePtr->applyCollisionData(it->second); + } + + if (!ratePtr->hasCrossSectionData()) { + throw CanteraError("addCollision", + "ElectronCollisionPlasmaRate requires either inline cross-section data " + "or a valid 'collision' reference."); + } + + if (!ratePtr->target().empty() && ratePtr->target() != target) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' targets '{}', but reaction '{}' uses target '{}'.", + ratePtr->collisionName(), ratePtr->target(), + collision->equation(), target); + } + + string kindFromReaction = inferElectronCollisionKind(collision); + string kindFromCollision = ratePtr->kind(); + + bool compatibleKind = kindFromReaction == kindFromCollision; + + if ((kindFromReaction == "elastic" || kindFromReaction == "effective") && + (kindFromCollision == "elastic" || kindFromCollision == "effective")) { + compatibleKind = true; + } + + // Allow collapsed inelastic channels. + // Example: + // reaction: Electron + N2 => Electron + N2 + // collision: kind: excitation, product: N2(rot) + if (!compatibleKind && + (kindFromReaction == "effective" || kindFromReaction == "elastic") && + kindFromCollision == "excitation") { + compatibleKind = true; + } + + if (!compatibleKind) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' has kind '{}', but reaction '{}' is inferred as '{}'.", + ratePtr->collisionName(), kindFromCollision, + collision->equation(), kindFromReaction); + } + m_collisions.emplace_back(collision); m_collisionRates.emplace_back( std::dynamic_pointer_cast(collision->rate())); @@ -457,7 +860,6 @@ void PlasmaPhase::addCollision(shared_ptr collision) updateInterpolatedCrossSection(i); // Set up data used by Boltzmann solver - auto& rate = *m_collisionRates.back(); string kind = m_collisionRates.back()->kind(); if ((kind == "effective" || kind == "elastic")) { @@ -475,11 +877,50 @@ void PlasmaPhase::addCollision(shared_ptr collision) m_kInelastic.push_back(i); } - auto levels = rate.energyLevels(); + + auto levels = ratePtr->energyLevels(); m_energyLevels.emplace_back(levels.begin(), levels.end()); - auto sections = rate.crossSections(); + auto sections = ratePtr->crossSections(); m_crossSections.emplace_back(sections.begin(), sections.end()); m_eedfSolver->setGridCache(); + +} + +// Infer the collision kind from the electron balance: producing electrons +// indicates ionization, consuming electrons indicates attachment, unchanged +// reactants/products indicate elastic scattering, and the remaining case is +// treated as excitation. +string PlasmaPhase::inferElectronCollisionKind( + const shared_ptr& collision) const +{ + const string eName = electronSpeciesName(); + + double nReactantElectrons = 0.0; + double nProductElectrons = 0.0; + + auto reactantElectron = collision->reactants.find(eName); + if (reactantElectron != collision->reactants.end()) { + nReactantElectrons = reactantElectron->second; + } + + auto productElectron = collision->products.find(eName); + if (productElectron != collision->products.end()) { + nProductElectrons = productElectron->second; + } + + if (nProductElectrons > nReactantElectrons) { + return "ionization"; + } + + if (nProductElectrons < nReactantElectrons) { + return "attachment"; + } + + if (collision->reactants == collision->products) { + return "elastic"; + } + + return "excitation"; } bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) @@ -526,7 +967,7 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() static const int cacheId = m_cache.getId(); CachedScalar last_stateNum = m_cache.getScalar(cacheId); - // combine the distribution and energy level number + // combine the distribution and energy level number to get a single cache variable. int stateNum = m_distNum + m_levelNum; vector interpChanged(m_collisions.size()); @@ -553,7 +994,17 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) { - // @todo exclude attachment collisions + if (m_elasticElectronEnergyLossCoefficients.size() != nCollisions()) { + m_elasticElectronEnergyLossCoefficients.resize(nCollisions(), 0.0); + } + + const string kind = m_collisionRates[i]->kind(); + + if (kind == "attachment") { + m_elasticElectronEnergyLossCoefficients[i] = 0.0; + return; + } + size_t k = m_targetSpeciesIndices[i]; // Map cross sections to Eigen::ArrayXd @@ -587,6 +1038,12 @@ double PlasmaPhase::elasticPowerLoss() // collisions (inelastic recoil effects). double rate = 0.0; for (size_t i = 0; i < nCollisions(); i++) { + const string kind = m_collisionRates[i]->kind(); + + if (kind == "attachment") { + continue; + } + rate += concentration(m_targetSpeciesIndices[i]) * m_elasticElectronEnergyLossCoefficients[i]; } @@ -777,7 +1234,8 @@ double PlasmaPhase::jouleHeatingPower() const return sigma * E * E; // W/m^3 } -double PlasmaPhase::intrinsicHeating() +// actually be careful, this is inelastic + growth +double PlasmaPhase::inelasticPower() { // Joule heating: sigma * E^2 [W/m^3] const double qJ = jouleHeatingPower(); @@ -787,8 +1245,92 @@ double PlasmaPhase::intrinsicHeating() double qElastic = elasticPowerLoss(); checkFinite(qElastic); - return qJ + qElastic; + return qJ - qElastic; +} + +double PlasmaPhase::intrinsicHeating() +{ + // Joule heating: sigma * E^2 [W/m^3] + const double qJ = jouleHeatingPower(); + checkFinite(qJ); + + checkVibrationalReservoirMoleFractions(); + + return qJ; +} + +void PlasmaPhase::updateVibrationalReservoirSpecies() +{ + m_vibrationalReservoirSpecies.clear(); + + for (size_t k = 0; k < nSpecies(); k++) { + string baseName; + const string& reservoirName = speciesName(k); + + if (!isVibrationalReservoirName(reservoirName, baseName)) { + continue; + } + + size_t kBase = speciesIndex(baseName, false); + + if (kBase == npos) { + warn_user("PlasmaPhase::updateVibrationalReservoirSpecies", + "Species '{}' matches the fictive vibrational-reservoir naming " + "convention, but the corresponding base species '{}' was not found. " + "This species will not be monitored as a vibrational reservoir.", + reservoirName, baseName); + continue; + } + + VibrationalReservoirSpecies reservoir; + reservoir.reservoirIndex = k; + reservoir.baseSpeciesIndex = kBase; + + m_vibrationalReservoirSpecies.push_back(reservoir); + } + + m_vibrationalReservoirSpeciesNeedUpdate = false; } +void PlasmaPhase::checkVibrationalReservoirMoleFractions() +{ + if (m_vibrationalReservoirSpeciesNeedUpdate) { + updateVibrationalReservoirSpecies(); + } + + for (const auto& reservoir : m_vibrationalReservoirSpecies) { + const size_t kReservoir = reservoir.reservoirIndex; + const size_t kBase = reservoir.baseSpeciesIndex; + + const double Xv = moleFraction(kReservoir); + const double Xb = moleFraction(kBase); + + const double pool = Xv + Xb; + + if (pool <= m_vibrationalAbsoluteMoleFractionThreshold) { + continue; + } + + const double reservoirFraction = Xv / pool; + + if (reservoirFraction > m_vibrationalMoleFractionThreshold) { + const string& reservoirName = speciesName(kReservoir); + const string& baseName = speciesName(kBase); + + writelog("Warning: fictive vibrational reservoir species '{}' contains " + "{:.3e} of the total '{}' pool. " + "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " + "Chemistry involving '{}' may be affected because part of the " + "material is stored in an inert vibrational reservoir.\n", + reservoirName, + reservoirFraction, + baseName, + reservoirName, Xv, + baseName, Xb, + m_vibrationalMoleFractionThreshold, + baseName); + } + } +} } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index b4d120099db..332ab9cb111 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -14,6 +14,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index fbe699b3889..d03da101282 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -13,6 +13,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h"